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INTRODUCTION 


One  of  the  difficulties  faced  in  the  analysis  of  structural/geo¬ 
technical  systems  is  numerically  modeling  a  semi - inf inite  domain.  The 
finite  element  method  (FEM)  provides  an  effective  numerical  tool  for 
modeling  complex  constitutive  relations  which  characterize  soil  behaviors 
however,  since  the  FEM  is  domain-based,  it  has  difficulty  modeling  the 
semi-infinite  domain.  Simply  truncating  the  domain  at  some  arbitrary 
distance  can  result  in  the  reflection  of  strain  energy  from  the  artifi¬ 
cial  boundary  which  should  continue  to  be  transmitted  into  the  semi¬ 
infinite  domain.  The  reflected  energy  can  contaminate  the  near  field 
region  where  the  solution  is  of  interest. 

Considerable  effort  has  been  devoted  to  developing  "silent"  or 
transmitting  boundaries  in  the  frequency  domain  which  allow  for  the 
outward  radiation  of  energy  (see  e.g.,  Waas,  1972;  Lysmer  and  Waas, 

1972;  Kausel,  Roesset,  and  Wass,  1975;  Bettess  and  Zienkiewicz,  1977; 
Karasudhi  and  Rajapakse,  1984).  However,  frequency  domain  solutions  are 
limited  to  linear  analyses. 

Lysmer  and  Kuhlemeyer  (1969)  developed  viscous  boundaries  in  an 
early  attempt  to  provide  a  "silent  boundary."  These  boundaries  are 
simple  and  can  be  used  in  both  time  and  frequency  domains.  However, 
they  act  as  perfect  absorbers  only  for  a  limited  class  of  problems. 

Cohen  and  Jennings  (1983)  give  an  excellent  review  of  the  different 
efforts  made  to  provide  a  "silent  boundary." 

In  recent  years,  the  boundary  element  method  (BEM)  has  become  a 
strong  candidate  for  use  in  the  analysis  of  structural/geotechnical 
systems  since  it  inherently  natisf ie.->  *■!>»•  r»«d'fntion  boundary  condition. 

In  this  study  we  concentrate  on  coupled  solution  approaches  where  boun¬ 
dary  element  methods  are  combined  with  the  finite  element  method.  The 
boundary  element  method  is  used  to  model  the  semi  -  inf inite  or  infinite 
domain.  In  particular,  we  consider  BEM  formulations  based  upon  the 
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Stokes  solution  --  the  analytical  solution  for  a  time  varying  con¬ 
centrated  load  in  an  infinite  domain. 

Objectives 

The  objectives  of  this  report  are  to  provide:  (1)  the  necessary 
analytical  and  numerical  background  for  time  domain  boundary  element 
methods  based  on  Stokes's  solution,  and  (2)  a  coupling  algorithm  for 
combining  the  BEM  with  the  FEM. 

The  ultimate  objective  of  this  research  is  to  determine  an  effec¬ 
tive  numerical  scheme  for  modeling  nonl inear,  dynamic  structural/geo¬ 
technical  problems  which  have  infinite  domains. 

Background 

The  use  of  integral  equation  formulations  in  the  analysis  of  tran¬ 
sient  phenomena  in  solids  and  fluids  dates  back  over  100  years  to  the 
Helmholtz-Kirchof f  integral  formula,  according  to  Manolis  (1984).  This 
formula  is  the  mathematical  description  of  Hygens'  principle  (Baker  and 
Copson,  1939). 

Though  integral  equation  statements  of  wave  propagation  phenomena 
have  existed  for  many  years,  their  numerical  approximation  has  occurred 
relatively  recently.  Among  the  early  efforts  were  the  works  of  Friedman 
and  Shaw  (1962)  and  Chen  and  Schweikert  (1963)  in  acoustics,  and  Banaugh 
and  Goldsmith  (1963)  in  steady-state  el astodynamics .  Cruse  and  Rizzo 
(1968)  and  Cruse  (1968)  were  the  first  to  apply  the  BEM  to  transient 
elastodynamic  problems.  Their  papers  considered  the  BEM  in  conjunction 
with  a  Laplace  transformation  to  solve  a  half-plane  wave  propagation 
problem.  Other  researchers  considered  methods  using  the  Fourier  trans¬ 
formation  instead  of  the  Laplace  transformation  (see  e.g.,  Niwa,  Kobayashi, 
and  Azuma,  1975;  Niwa,  Kobayashi,  and  Fukui,  1976).  Shaw  (1985)  gives 
an  overview  of  many  different  BEM  formulations  in  elastodynamics . 

The  transformation  methods  have  been  highly  developed,  but  to 
address  problems  where  the  FEM  region  is  characterized  by  nonlinear 


2 


behavior  we  must  consider  time  domain  approaches.  Three  time  domain 
approaches  have  been  pursued.  They  differ  in  one  of  the  most  basic 
aspects  of  the  BEM  --  the  fundamental  singular  solution.  The  two  sim¬ 
pler  approaches  are  based  on  fundamental  singular  solutions  for  elliptic 
partial  differential  equations  instead  of  the  actual  hyperbolic  equations. 

One  of  the  approaches  approximates  the  time  derivatives  by  finite 
difference  (see  e.g.,  Brebbia  and  Walker,  1980).  The  resulting  partial 
differential  equation  is  solved  by  the  boundary  element  method  at  each 
time  step.  This  method  requires  domain  integrations  of  previous  time 
step  displacement  fields  because  of  the.  finite  difference  approximation; 
thus  it  would  appear  to  have  limited  application  for  problems  with  semi- 
infinite  or  infinite  domains. 

Another  approach  uses  the  Kelvin  solution  of  elastostatic  BE  formu¬ 
lations.  The  main  idea  behind  this  approach  is  to  approximate  the  inertia 
effects  by  expanding  the  displacement  field  throughout  the  domain  in 
terms  of  a  special  set  of  expansion  functions.  This  approximation  is 
only  used  for  the  inertia  terms  and  the  expansion  functions  are  neces¬ 
sarily  simple  allowing  the  domain  integration  to  be  written  in  terms  of 
boundary  integrations.  Nardini  and  Brebbia  (1982)  first  developed  this 
approach  and  then  followed  with  several  applications  (Nardini  and  Brebbia, 
1983  and  1986). 

Later  Ahmad  and  Banerjee  (1986)  used  the  concepts  of  complementary 
functions  and  particular  integrals  to  solve  the  free  vibration  problem 
in  elastodynamics .  Though  the  derivation  of  their  approach  is  different 
than  the  previous  work  by  Nardini  and  Brebbia,  the  resulting  system  of 
equations  is  almost  identical.  The  reader  is  advised  to  verify  the  com¬ 
parisons  made  in  this  later  work.  Both  of  these  approaches  have  been 
applied  effectively  to  bounded  domain  problems;  however,  I  am  unaware  of 
any  work  which  has  successfully  applied  them  to  unbounded  domain  problems. 

The  most  rigorous  time  domain  approach  is  based  upon  the  fundamental 
solution  of  elastodynamics  --  the  Stokes  solution.  This  appears  to  be 
the  best  suited  approach  for  handling  infinite  domains.  The  fundamental 
solution  satisfies  the  radiation  boundary  conditions,  and  for  quiescent 
initial  conditions  with  no  body  forces,  the  method  requires  no  domain 
integrations.  However,  its  formulation  includes  a  convolution  of  the 
time  variable  and  thus  appears  to  be  computationally  intense. 


Attempting  to  provide  a  more  efficient  solution  to  problems  with 
infinite  domains,  researchers  (see  e.g.,  Geers,  1983)  have  developed 
simplified  BEN  formulations  known  as  doubly  asymptotic  approximations 
(DAA).  Mita  and  Luco  (1987)  give  an  overview  of  the  different  BEM  and 
coupled  FEM-BEM  approaches  which  have  been  used  to  study  the  dynamics  of 
embedded  foundat ions . 

Scope 


The  scope  of  this  work  is  limited  to  a  discussion  of  time  domain 
BEMs  which  are  based  on  the  use  of  the  Stokes  solution.  Both  analytical 
and  numerical  aspects  of  the  formulations  are  presented.  The  analytical 
roots  of  the  formulations  are  presented  in  some  detail;  this  provides  a 
strong  background  for  future  implementation.  The  emphasis  is  on  the 
direct  BEM  formulation  but  the  indirect  BEM  is  also  presented.  The 
coupling  algorithm  treats  the  BEM  region  as  a  nonlinear  boundary  con¬ 
dition  for  the  FEM  region.  Equilibrium  and  displacement  compatability 
on  the  interface  between  the  regions  is  satisfied  in  a  nodal  sense  using 
an  iterative  scheme. 


THEORY 

This  section  presents  the  analytical  basis  for  time  domain  BEMs 

based  on  the  Stokes  solution.  The  Nomenclature  at  the  end  of  this  report 

presents  an  additional  explanation  of  some  of  the  symbols  used  in  this 

section.  For  additional  detail  the  interested  reader  could  see,  for 

example,  Eringen  and  Suhubi's  (1974,  1978)  excellent  texts  on  elasto- 

dynamics.  The  later  text  is  the  main  reference  for  this  section. 

The  motivation  for  coupling  the  FE  and  BF,  methods  to  solve  a 

boundary-initial  value  problem,  is  to  apply  each  method  to  that  portion 

of  the  domain  for  which  it  is  best  suited  fsee  Figure  1).  The  finite 

F 

element  method  in  this  formulation  will  be  applied  to  a  subdomain  ft 
which  includes  all  portions  of  the  domain  that  are  nonlinear,  inhomo¬ 
geneous,  or  anisotropic.  The  boundary  element  method  in  this  lorrtiu- 
lation  will  be  applied  to  a  subdomain  ft^  which  will  be.  idealized  as 
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consisting  of  a  homogeneous,  isotropic,  linear  elastic  material  subjected 
to  small  displacements  and  deformations.  The  strength  of  the  BEM  is  its 
inherent  ability  to  model  infinite  domains.  The  following  section  presents 

g 

the  governing  equations  for  ft  --  the  equations  of  linear  elastodynamics . 

g 

For  brevity,  ft  will  be  denoted  as  ft  since  the  following  sections  only 
address  the  boundary  element  subdomain. 


Elastodynamics  Problem 


Field  Equations.  The  governing  equations  of  motion  are  given  by: 


0 .  .  .  +  <l> .  =  pu . 

Ji.J  i  i 


a  .  .  =  o  .  . 

ij  Ji 

and  the  compatibility  or  kinematic  equations  are.  given  by: 


£,,  =  r  U,  ,  ,  U.  ,  (Z) 

ij  2  \  i,j  + 

The  BE  subdomain  is  assumed  to  be  an  isotropic  homogeneous  linear 
elastic  material.  Thus  the  governing  constitutive  relations,  generalized 
Hook's  law,  are  given  as: 

a .  ,  =  2jje  .  .  +  Xe  6 .  .  (3) 

ij  ij  kk  ij 

where  X  and  y  are  Lame's  constants  expressed  in  terms  of  Young's  modulus 
(E)  and  Poisson's  ratio  (v)  as 


Ev 

(1+v)  ( l-2v)  ’ 


E_ 

2(j+v) 


Alternatively,  the  equations  of  motion  (Equation  1),  compatibility 
relations  (Equation  2),  and  constitutive  relations  (Equation  3)  can  be 
combined  into  a  linear  system  of  hyperbolic,  differential  equations  in 
terms  of  displacement  as: 


uui.jj  +  <x+u)uj,ji  =  piii 


These  are  the  well  known  Navier  equat  ions  of  elastodynamics .  These 
equations  are  often  written  in  a  slightly  different  form  which  introduces 
some  of  the  complications  that  occur  with  vector  hyperbolic  equations. 

The  Stokes -He  Imho  1 tz  resolution  theorem  states  that  every  sufficiently 
smooth  vector  field  f(x,t)  may  be  decomposed  into  irrotational  and 
solenoidal  parts;  that  is,  it  admits  the  representation: 

f  =  Vf  +  V  x  F  (6) 

where  the  first  term  is  curl  free  and  the  second  term  is  divergence  free. 
Applying  this  to  both  the  displacement  and  body  force  vector  fields  we 
can  write: 


P 

u 


V  0  , 


S 


u 


V  x  X 


and 

T  =  Vf  +  V  x  F 


(7a) 


(7b) 


where  the  scalar  valued  function,  0,  and  the  vector  valued  function,  X, 

are  the  Lame  potentials  (F.ringin  and  Suhubi,  1975).  Substituting  these 

relationships  into  Equation  5  can  show  that,  the  Navier  equations  are 
P  S 

satisfied  if  u  and  u  satisfy: 


?op  p  c 

CpV  u  +  Vf  =  ii  ,  (TV  u  +  V  x  F  =  ii  (8) 

respectively.  In  the  absence  of  body  forc.es,  these  equations  yield  the 
familiar  vector  wave  equations  for  irrotational  and  equivoluminal  wave 
propagation,  respectively.  These  waves  propagate  with  velocities  given 
by: 


X+2y 


\J  P 


(9) 


Noting  that  X  and  y  are  positive  implies  that  C  >  C  and,  thus, 

i  u 

the  rationale  for  the  subscripts:  P  for  'primary'  corresponding  to  the 
faster  wave,  and  S  for  'secondary'  corresponding  to  the  slower  wave.  In 
addition  to  irrotational  and  primary,  the  P  waves  are  referred  to  as: 
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dilational,  pressure,  compress iona  l ,  and  longitudinal  waves.  In  addition 
to  equivoluminal  and  secondary,  the  S  waves  are  referred  to  as:  distor- 
tional,  shear,  and  transverse  waves. 

Navier's  equations  are  otten  written  in  terms  of  these  wave  speeds 
as  : 


cln. 

S  i. 


JJ 


K  •  cs2)' 


j.  j1 


r  = 

i 


no) 


The  equations  of  motion,  kinematic  relations,  and  constitutive  rela¬ 
tions  govern  the  problem  ov  r  the  domain  except,  for  singular  surfaces 
where  special  jump  conditions  must  be  satisfied.  The  wave  front  of  a 
shock  wave  is  an  example  of  a  singular  surface.  For  additional  details 
on  singular  surfaces  and  the  corresponding  jump  conditions  see  Eringen 
and  Suhubi  (1974). 


Initial  Conditions.  The  initial  conditions  give  the  displacement 
and  velocity  field  throughout  the  domain  ns: 


ui(x,0) 

x  e  ft 

(11a) 

ui(x, 0) 

=  voI(x) 

X  £  fi 

(lib) 

Boundary  Conditions. 

The  boundary  conditions  are  given  by: 

ui(x, t) 

=  ui(x,t) 

X  E  ,  t  E  |  0,») 

(12a) 

ti(x, t) 

=  E^x.t) 

X  E  fT  ,  t.  E  [  0 , » ) 

(12b) 

where  u(x,t)  and  £(x,t)  are  prescribed  distributions  of  boundary  displace¬ 
ments  and  tractions  as  a  function  of  time.  The  simple  notation  does  not 
imply  that  the  boundary  conditions  are  mutually  exclusive;  the  fully-mixed 
boundary  value  problem  is  addressed. 

Radiation  Boundary  Conditions.  if  the  domain  is  unbounded,  physics 
places  constraints  on  the  behavior  of  the  fields  at  infinity.  Physical 
reasoning  suggests  that  if  the  applied  loading  is  restricted  to  a  finite 
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region,  waves  propagating  from  infinity  into  the  interior  of  the  domain 
should  not  exist. 

Radiation  boundary  conditions  are  attributed  to  Sommerfeld  (1949). 

He  suggests  the  necessity  of  such  a  condition  in  the  following  discussion: 

"...  oscillation  problems  (in  constrast  to  potential  problems)  are 
not  determined  uniquely  by  their  prescribed  sources  in  the  finite 
domain.  This  paradoxical  result  shows  that  the  condition  of 
vanishing  at  infinity  is  not  sufficient,  and  that  we  have  to 
replace  it  by  a  stronger  condition  at  infinity.  We  call  it  the 
condition  of  radiation:  the  sources  must  be  sources,  not  sinks  of 
energy.  The  energy  which  is  radiated  from  the  sources  must  scatter 
to  infinity;  no  energy  may  be  radiated  from  infinity  into  the  pre¬ 
scribed  singularities  of  the  field." 


Regardless  of  whether  we  are  seeking  a  closed  form  or  a  numerical 
solution,  the.  radiatior.  condition  provides  the  same  essential  property 
to  our  solution  --  uniqueness. 

It  can  be  proven  (Eringen  and  Suhubi ,  1975)  that  the  radiation  con¬ 
ditions  of  elastodynamics  are  direct  consequences  of  the  radiation  con¬ 
ditions  on  wave  equations  (Equation  8).  The  radiation  conditions  of 
elastodynamics  can  be  stated  as: 


lim  r 


0 


(13a) 


lim  r 
r+o© 


0 


(13b) 


where  t  and  t'  are  traction  vectors  on  a  sphere  of  radius  r,  due  to  the 

P  S 

displacement  components  u  and  ul  ,  respect  ively.  These  conditions  are 
sufficient  to  guarantee  that  at  infinity  there  will  only  be  an  outward 
flow  of  energy;  that  is,  reflections  are  eliminated. 
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Integral  Equation  Formulations 


This  section  presents  the  integral  equations  on  which  both  the  direct 
and  indirect  BF.Ms  are  based.  Only  the  equations  for  the  direct  formula¬ 
tion  are  developed  in  detail.  The  development  of  the  direct  formulation 
in  elastodynamics  strongly  parallels  the  development  in  elastostatics ; 
the  two  major  components  are  a  fundamental  singular  solution  and  a  re¬ 
ciprocal  theorem.  Kelvin's  solution  is  replaced  by  Stokes's  solution, 
and  Betti's  reciprocal  theorem  is  extended  to  Graffi's  theorem. 

Elastodynamic  State.  Prior  to  presenting  Stokes's  solution  and 
deriving  Graffi's  theorem,  it  is  useful  to  have  the  definition  of 
elastodynamic  state  as  given  by  Eringen  and  Suhubi  (1975,  Section  5.7). 

Let  R  be  a  spacial  region  with  boundary  T,  and  T  a  time  interval. 

If  u  and  a  are,  respectively,  a  vector-valued  and  a  symmetric  second- 
order  tensor-valued  function  defined  on  52  x  T,  we  call  the  ordered  pairs 
y  =  [u,d]  an  elastodynamic  state  on  52  x  T  with  the  displacement  field  u 
and  stress  field  a,  corresponding  to  a  body  force  density  Y,  mass  density 
P,  irrotational  wave  speed  Cp,  and  equivoluminal  wave  speed  Cg,  provided 
that : 

(a)  u  e  C2,2(S2  x  T),  u  e  C1,1(r  x  T) ,  o  e  C°’°(52  x  T) , 

as  I  E  C°’°(52  x  T),  p  >  0,  Cp  >  2//3  Cg  >  0  (14) 

(b)  u,  O,  Y,  p,  Cp,  ar.d  C  satisfy  the  governing  Equations  1,  2, 

and  3.  b 

The  class  of  all  elastodynamic  states  satisfying  the  above  conditions 
is  denoted  by  E  where  we  write: 

y  t  E  (f,p,Cp,Cs;  S2  x  T)  (15) 

When 

T  =  T  and  u  =  0  on  52  x  T 

we  refer  to  y  as  an  elastodynamic  state  of  quiescent  past  and  write: 
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(16) 


y  t  Eo(T,Plcp,cs;  n) 

See  Nomenclature  at  end  of  report  for  a  description  of  the  notation 
used  in  the  continuity  conditions  of  (a)  above. 

Stokes's  Solution.  Stokes's  problem  consists  of  an  infinite  domain 
subjected  to  a  concentrated  load  at  a  point  £  which  is  fixed  in  direction 
e.  but  has  an  arbitrary  time  varying  magnitude.  That  is,  we  seek  the 
solution  in  an  infinite  domain  where  the  body  force  is  given  as: 

^(x.t)  =  f(t)  6(x-£)  e^  (17) 

The  solution  is  the  fundamental  singular  solution  of  elastodynamics 
and  is  originally  due  to  Stokes  (1849)  (also  see  Love,  1944,  Section  212 
or  Eringen  and  Suhubi,  1975,  Section  5.10).  The  displacement  at  a  position 
x  and  time  t  is  given  as: 

u . (x, t)  =  u . 

i  —  i 

and 

uij(x,t;^| f) 

u^(x,  t;£|  f )  = 


j  ( x » t ;  £  I  f )  e. 


(18a) 


u*j\(x,  t;£|  f)  +  u ^  ,(x,  t;£  |  f ) 


(18b) 


1 

4trpr 


,  3r . r  . 


.  )  j  *  f(t-Xr)  dX 


+  f't  _  e  ) 

rV  '  Cp) 


(18c) 
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-1 


uij(x,t;i;j  f) 


4Trpr 


3r  .  r  . 

LA 


ij 


X  f(t-Xr)  dX 


r .  r  . 

~ti  f|t 

r  cs 


<■  -  y  *  %  -(•  ■  y 


(18d) 


r.  =  x.  -  l. 

1  l  i 


r  =  r  .  r  . 

l  l 


( 18e) 
(18f) 


P  S 

u_  and  u  are  the  irrotational  and  equivoluminal  parts  of  the  Stokes 
tensor  u_.  Substituting  this  result  into  the  kinematic  (Equation  2) 
and  constitutive  equations  (Equation  3)  one  can  obtain  the  stress  at 
(x,t)  as  (Eringen  and  Suhubi,  1975): 


OjjCx.t)  =  o . jk(x,t;£| f)  eR 


(19a) 


where 


f)  =  a^k(x,t;4|f)  +  o^.k(x,t;5|  f) 


(19b) 


f) 


1 

' 

5r.r.r,  6  .  .  r,  +  6,,r.  +  6,,r, 

6  cs 

i  ,i  k  i.i  k  ik  j  jk  i 

4irr2 

3  r 

r 

-1 


P  2 

r  2Cs 

X  f(t-Xr)  dX  -  ~ 

6r.r.r,  6.  .r.  +  5,.  r.  +6..  r. 

l  j  k  ij  k  ik  j  jk  i 

3  r 

0  ^P 

r 

-  k)  -  (‘  - 2  J)  N‘  •  k) +  k  *('  -  k) 


2r  ,  r  ,r,  C 


Ol*5  .  I_) 

r2  C3  ' 


(19c) 
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4irr 


6  C 


5r  ,  r  .  r, 
.-1  J  k_ 


6  .  .  r.  +  6  ..  r .  +  6  ..  r  , 
ij_  k  lk  j  jk  l 


s 

j  X  f(t-Xr)  dX  +  2 

6r.r,r,  6,.r.  +  6.,r.  +  6.,r. 

i  J  k  ij  k  lk  j  jk  l 

3  r 

0 

!  r  J 

(t  r  I  SikrJ  +  ‘jk1-! 

.  JL.)  +  JL  .  £_V 

r  cs> 

T  cs)  +  csf(t  Cgjj 

+ 


fvjf* 

r2C„ 


(19d) 


Note  that  with  the  above  notation  we  can  consider  u..  and  a...  as 
operators  on  the  time  variation  of  the  concentrated  load  f(t).  The 
quantities  of  the  form  (t-r/C)  are  referred  to  as  retarded  times;  the 
effect  of  the  temporal  variation  of  the  load  is  retarded  by  the  time  it 
takes  a  wave  to  travel  from  the  source  to  the  field  point. 

The  elastodynamic  state,  which  corresponds  to  Stokes's  solution 
with  the  concentrated  force  acting  parallel  to  the  x^-axis,  is  referred 
to  as  the  Stokes  state  and  denoted  (Eringen  and  Suhubi,  1975,  Section 
5. 10)  by: 


t/k(x>t;£|f)  =  |yk(x,t;£|f),  ok(x,t;£|f)J 

where  the  vectors  u^  and  second-order  tensors  ( k= 1,2,3)  are: 

%  -  [uik]-  °k  -  hjk] 


(20a) 


(20b) 


In  particular,  we  will  often  need  to  refer  to  a  class  of  Stokes's 
states  that  have  a  quiescent  past.  Consider  the  following  definition 
given  by  Wheeler  and  Sternberg  (1968): 
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Let  £  £  and  f  be  a  twice  continuously  differentiable  function  of 
time  that  vanishes  on  T  ,  and  p,  Cp,  Cg,  satisfy  Equation  14.  We  call 
the  state  t/^Cx. t: ;£|  f )  =  [u^,o^]  defined  on  x  T  by  Equations  20,  18, 
and  19  the  Stokes  state  of  quiescent  past  for  a  concentrated  force  at  £ 
parallel  to  the  x^-axis  corresponding  to  the  force  function  f(t)  and  to 
the  material  constants  p,  C  ,  C  . 

r  o 

If  we  assume  Equation  18  is  valid  for  an  impulse  force,  i.e., 
f(t)  =  6(t-t) ,  the  substitution  leads  to  the  free  space  Green's 
function: 


Gij(x.t;£.'0  =  u  [x,t;£| 6(t-t) ] 

.  th 


(21a) 


So  this  expression  gives  the  i  component  of  the  displacement  field 
at  (x,t)  due  to  the  component  of  a  concentrated  impulse  acting  at 
(£,t).  In  some  contexts  the  Green's  function  is  written  for  t  =0  as: 


(21b) 


By  considering  the  "sifting  property"  of  the  6-function,  integration 
of  Stokes's  solution  gives  the  Green's  function  as: 


Gij<-’t;£,T) 


where  t '  =  t  -  t 


t' 


4irpr 


](~rj  -  Tj)  -(*•  - 1)  - »(*'  •  k 


6. 


r .  r . 

+ 


l  6(t*  ‘  cj  ■  Cc  6(t’  "  cj 


(22) 
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The  above  result  for  the  displacement  field,  when  substituted  into 
the  kinematic  and  constitutive  relations,  gives  the  stress  field  as: 


=  T, jk(x,t;£,t)  ek(£) 


(23a) 


where 


Tijk(-’t;£,T) 


1  1 

-  6  C> 

r.r.r,  6.  .r,  +  6.,  r.  +  6.,  r. 

r  1  J  k  1.1  k  lk  j  jk  l 

47rr2 

6  2 
r 

k>  0  - 

3  r 

r 

H(c'  - 1)  -  H(c’  -  k 


+  2 


r.r.r,  6 .  .r,  +  6.,r.+6.,r. 

6  -1  ..J_k  _  _U  k  _ jiL  1 


+  2 


(*'  5(‘’  -  §;) 

*(*'  '  t)  4-  »(f  - 1) 


r  .  r .  r, 

.  1  J  k 


r2  C, 


9M1  -  *?-)[•(*■- *{*'- U 


6  ..  r  .  6  ,  .  r 
ik  j  ij  k 


-szl  +  M1'  -k) 


(23b) 


By  Cauchy's  formula  one  can  obtain  the  relationship  for  the  traction 


on  a  plane  with  unit  normal  n,  as: 

i 


ti(x,t)  =  Fik(x,t;£,i)  ek(£) 


(24a) 
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where 


Flk(x,t;£,i) 


T1jk(x»t;i.'r)  n  (x) 


(24b) 


It  should  be  noted  that  many  equivalent  forms  of  Equations  22  and 
23  are  given  in  the  literature.  The  above  classical  solutions  to  the 
equations  of  elastodynamics  are  an  essential  element  to  the  time  domain 
boundary  element  methods  of  this  study.  The  second  essential  element  is 
the  dynamic  reciprocal  theorem  considered  in  the  following  section. 

Dynamic  Reciprocal  Theorem.  In  elastostatics  the  Betti-Rayleigh 
reciprocal  work  theorem  provides  a  relationship  between  two  distinct 
equilibrium  states.  The  dynamic  reciprocal  theorem  can  be  written  by 
applying  the  D'Alembert  principle  --  including  the  inertia  terms  as  a 
part  of  the  body  force.  By  this  approach,  Rayleigh  (1873)  obtained  a 
theorem  (Love,  1944,  Section  121)  from  which  Graffi's  theorem  (1946-1947) 
can  be  deduced;  we  will  take  this  approach  in  the  following  discussion. 

For  another  proof  of  Graffi's  theorem  see,  for  example,  Eringen  and  Suhubi 
(1975,  Section  5.8). 

The  following  proof  of  Graffi's  theorem  is  included  for  more  than 
the  sake  of  completeness.  A  numerical  method  cannot  be  properly  used, 
let  alone  be  implemented,  without  a  basic  understanding  of  the  underlying 
analysis.  The  development  which  follows  also  uses  notation  common  to 
much  of  the  literature. 

Consider  two  distinct  elastodynamic  states: 


V  =  [u,o]  e  E(t, p ,Cp, Cg ;  8  x  T+)  (25a) 

y*  =  [u*,d*]  e  E(T*,p,Cp,Cs;  8  x  T+)  (25b) 

defined  on  8  with  initial  conditions: 

u(x,0)  =  u  (x)  u(x,0)  =  v  (x)  (26a) 

—  —  — o  —  —  —  o 

in  8 

u  (x,0)  =  uq(x)  M  (x,0)  =  yQ(x)  (26b) 


The  elastodynamic  equilibrium  equations  for  each  state  are  given  as: 
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0 


(27a) 


o  ,  .  .  +  $  ,  — 

ij ,  J  i 


*  * 

a  .  .  .  +  e  ,  =  0 

iJ,J  i 


in  ft 


(27b) 


where 


32u, 


Pi 


* 


(V  l-i  ,  \ 

h--i) 


(28a) 


3x 


.2  * 
3  u . 


‘  +1  ' 


(28b) 


Betti's  reciprocal  theorem  for  two  distinct  equilibrium  states  is 
given  as  (Sokolnikof f ,  1956,  Section  109): 


f  *  (  *  {  *  f  * 

I  t.u.  dA+  I  3.u.  dV  =  I  t.u,  dA+  I  fl.u. 

Jii  Jit  J  ii  Jii 


dV 


(29) 


n 


Note  that  we  can  substitute  Equation  28  into  Equation  29.  Since 
this  expression  is  true  for  all  x  e  T+,  we  can  integrate  from  0  to  t  to 


obtain: 

t 


t  .2  * 

3  u . 


f  f  t*u.  M  dt  +  f  f  pY*U.  dV  dp  -  f  f  P  — i  », 

a  r  a  n  a 


dV  dx  = 


o  r 

t 


o  n 

t 


o  n 


t  .2 

3  u , 


I  j  d*  d,  ♦  J  f  PT  ,u*  dV  dx  -  |  |  P  - y-  u*  dV  dx  (30) 

or  on  o  n  3t 

Integration  by  parts  of  the  time  derivative  terms  gives: 
t  2  * 

r  3 

J  - -  u^(x, x )  dx  =  ui(x,t)  u.(x,t)  -  vQl(x)  uoi(x) 


-  J  u^(x,x)  u.(x,x)  dx 


Ola) 
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Vr  tV 

u.(x,t)  u.(x,t)  -  v  .(x)  u  .(x) 

1  —  1  —  Ol  —  Ol  — 

t 

* 

u.(x,x)  u.(x,x)  dx  (31b) 

0 

where  =  du^/dx .  Substituting  these  results  into  Equation  30  and  showing 
the  independent  variables  explicitly,  the  following  form  of  the  reciprocal 
theorem  of  elastodynamics  is  obtained: 


t  „2  ,  , 

r  3  u.(x,t)  * 

I  - 2 -  ui(x,r)  dx  = 


t 

I  f  ti(x,T)  ui(x,i)  dA(x) 

o  r 


t 

dt  +  j  j  p(x)  T^Cx,!)  u.(x,t)  dV(x)  dt 
0  Q 


-  I  p(x)  ui(x,t)  u1(x,t)  -  VQi(x)  uoi(?)  dV(x) 
Q 


t  t 

~  j  j  u*(x.T)  dA(x)  dt  +  j  j  p(x)  Ii(x,t)  ui(x,t)  dV(x) 

or  on 


dt 


-  J  p(x)  ui(x,t)  ui(x,t)  -  VQi(x)  uoi(x)  dV(x) 

n 


(32) 


Mansur  and  Brebbia  (1985)  derive  the  elastodynamic  equivalents  of 
Somigliana's  identities  from  this  relationship.  We  will  show  the  more 
classical  approach  obtaining  the  same  integral  equations  via  Graffi's 
theorem. 

In  obtaining  Equation  32,  we  tacitly  assumed  that  the  two  elasto¬ 
dynamic  states  occurred  at  the  same  time,  x.  Alternatively,  we  can,  in 
effect,  integrate  the  y*  state  from  t  down  to  0;  that  is,  the  reciprocal 
theorem  "compares"  the  y*  state  at  time  t'  =  t-x  with  the  y  state  at 
time  x.  Then  the  integration  by  parts  of  Equation  31  becomes: 
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t  2  * 

3  u.(x.t-x)  * 

- 7, - U.(x»T)  dx  =  U.(x,o)  Ii.(x,t)  -  u.(x,t)  u  .(x) 

/  1  1  1  —  1  “  Ol  _ 


/ 

0 


3x 


-  j  u^(x,t-x)  u^(x, x)  dx 


(33a) 


t  „2  , 

r  3  Ui(x,X) 

3t2 


k  *  k  k 

u.(x,t-x)  dx  =  u.(x,t)  u  .(x)  -  v  ,(x)  u.(x,t) 

1  ~  1  —  Ol  “  Ol  ~  1“ 


J  u.(x,x)  u.(x,t-x)  dx  (33b) 

0 


We  want  the 
to  be  taken  with 


time  derivatives  that  are  outside  of  the  time  integrations 
respect  to  t  instead  of  x.  Noting  that, 


k  k 

3u^(x,t-x)  3u^(x,  t-x) 

Tr  ~  =  it 

We  now  define: 

u.  =  3u./3t 

l  i 

Combining  the  results  of  Equation  33  with  Equation  30  and  rearranging 
terms  gives: 


j  j  t.(x.t-x) 

r  o 


u  (x,x)  dx 


t 

dA(x)  +  |  p(x)  | 
fl  0 


f  '(x,t-x) 


u.(x»t) 


dx 


k  k 

+  vQi(x)  ui(x,t)  +  uoi(x)  u  (x.t) 


dV(x) 
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t 

=  f  i  t^X.T) 
r  o 


* 

Ui(x, t-T) 


dx  dA(x) 


+  J  p(x) 

Q 


t 

J  ^(x.1)  u.(x,t-x)  dx 
0 


■fc  #  Vc 

+  v  .(x)  u.(x,t)  +  u  .(x)  u.(x,t) 


Ol 


Ol 


dV(x) 


(34) 


This  is  Graffi's  reciprocal  theorem  of  elastodynamics .  Thus  the 
dynamical  reciprocal  theorem  is  an  extension  of  Betti's  reciprocal  theorem 
of  elastostatics .  Most  authors  express  Equation  34  in  terms  of  Riemann 
convolutions  (see  Nomenclature)  as: 


/  dA(x)  +J  P(x) 

r  Q 


[OuiKx.t) 


*  Vr 

+  v  ,(x)  u.(x,t)  +  u  .(x)  u.(x,t) 

/"VI  —  1  — '  /"VI  —  1  —  ’ 


Ol 


Ol 


dV(x) 


J  (t^*ui](x,t)  dA(x)  +  J  p(x) 


[  ^ i*Ui 1 (-’ t) 


+  v  .(x)  u,(x,t)  +  u  .(x)  u.(x,t) 

Ol  _  1  —  Ol  -  1  ~ 


dV(x) 


(35) 


In  elastostatics,  an  integral  equation  statement  of  the  boundary 
value  problem  (Somigliana' s  identities)  is  derived  from  Betti's  recip¬ 
rocal  theorem  by  using  as  one  of  the  equilibrium  states  a  state  given  by 
the  fundamental  solution  (Kelvin's  solution).  In  the  following  section 
an  integral  equation  statement  of  the  elastodynamics  problem  is  derived 
from  Graffi's  reciprocal  theorem  using  Stokes's  solution  to  define  one 
of  the  elastodynamic  states. 


Integral  Equations  for  the  Direct  Boundary  Element  Method.  The 

integral  equation  statement  derived  in  this  section  is  an  exact  statement 
of  the  elastodynamics  problem.  The  numerical  approximation  of  this  state¬ 
ment  is  presented  in  a  later  section. 


19 


We  apply  Graffi's  theorem  to  two  e 1  as todynamic  states:  y  corresponds 

ic 

to  our  elastodynamics  problem  and  y  corresponds  to  the  Stokes  state  of 
quiescent  past  where  the  body  load  is  a  concentrated  impulse  load  applied 
at  t=0  (i.e.,  the  Green's  function  results  of  Equations  21  through  24). 
Graffi's  theorem  is  then  written  as: 


t 

I  J  F.k(x,t-T ;£)  ek(£)  u.(x,x)  dr  dA(x) 

r  o 


+ 


t 

f  f  fi(t-T)  6(x-C) 

n  o 


e  .  (£ )  ti  .  (  x ,  x  )  dx  dV(x) 

l  l 


t 

=||  t.(x,T) 

r  o 


6ik(x,t-x;g) 


ek(?)  dx  dA(x) 


t 

j  T.(x,x)  G.k(x,t-x;|)  ek(4)  dx 
0 


+  vo.(x)  G.k(x,t;5)  eR(4)  +  u^.fx)  G.k(x,t;£)  e_k(£) 


dV(x) 


(36) 


As  indicated  in  Equations  22  and  23,  the  times  t  and  x  always  occur 

in  the  Green's  function  (and  the  corresponding  higher  order  tensors)  as 

the  difference  t-x.  Physically,  the  response  of  the  domain  to  a  unit 

impulse  is  a  function  of  the  elaspsed  time  since  the  impulse  has  occurred. 

Mathematically,  the  Stokes  solution  (including  the  special  case  for  the 

Green's  function)  has  the  property  of  time  translation.  In  particular, 

let  B.,  represent  F.,  or  G.,  ,  then  we  have: 
lk  r  lk  ik 


Bik(x,t~T  ;£)  =  B,k(x,t-x  '.£.0)  =  B  jk(x,t;£,x)  (37) 

By  Equation  37  and  the  definition  of  the  Dirac  delta  "function," 
Equation  36  becomes: 
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where 


t 

I  J  |^ik(x,t;£,T)  t.(x,x)  -  Fik(x,t;4,T)  u^x.t)  di  dA(x) 
r  0 

t 

+  j  j  Gilc(x.t;^,T)  *.(x,t)  di  dV(x)  +  Ik(i,t; vq,Uo)  (38a) 

Q  0 


T^(4.t;vo,uo)  =  p  j  Jvq.(x)  G.k(x,t,^)  +  uq.(x)  G.k(x,t;£)  dV(x) 

Q 

Note  that  lP(£,t;v  ,u  )  accounts  for  the  effect  of  the  initial  con¬ 
ic  o  o 

ditions  on  the  displacement  field.  Two  alternative  forms  which  differ 
only  in  notation  are  given  below.  The  first  form  uses  the  convolution 
notation  and  is  often  written  in  the  literature  in  an  even  more  concise 
manner. 


uk(£,t)  =  /  {lGik*t.](x,t)  -  (F.k*u.)(x,t)  -  dA(x) 

r 

+  |  fGik**.J(x,t)  dV(x)  +  Ik(£,t;yo,uo)  (38b) 

The  second  form  is  expressed  in  terms  of  functional  operators  for 
«/k(x,  t;£)  f) ,  the  Stokes  state  of  quiescent  past  (F.ringin  and  Suhubi, 

1975,  Section  5.11).  It  is  given  as: 

=  j  |u°k(x. t .(x, t) J  -  t°k(x,t;£|u.(x,t)]  dA(x) 

r 

+  |  u°k(x.t;i|0i(x,t) |  dV(x)  +  I k ( 5 , t ; vq , uq )  (38c) 

Q 

where,  by  Cauchy's  formula: 

t°kIx,t;5|f]  =  o°jk[x,t;5| f]  n . (x) 
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The  equivalence  of  the  functional  operator  form  (Equation  38c)  to 
the  time  convolution  form  (Equation  38a  or  38b)  is  illustrated  by  Eringen 
and  Suhubi  (1975).  Equation  38  is  sometimes  referred  to  as  Love's  (1944) 
integral  identity  (Eringen  and  Suhubi,  1975,  give  additional  historical 
background).  Love  (1944)  considered  I^(£ ,  t  iv^u^)  in  more  detail  showing 
that  the  initial  displacement  and  velocity  at  a  point  x  affects  a  region 
bounded  by  two  spheres  traveling  at  wave  velocities  Cp  and  C  centered 

1  O 

at  x  (see  Eringen  and  Suhubi,  1975).  Equal  ion  38  is  the  elastodyr.amic 
counterpart  to  Somigliana's  first  identity  of  elastostatics  and  plays 
the  same  role  in  development  of  the  direct  BF,M  for  elastodynamics . 

The  equation  does  not  represent  a  solution  to  the  e lrs todynamic 
problem  since  the  traction  and  displacement  distribution  along  the 
boundary  are  only  partially  known.  (For  a  well  posed  problem,  only 
"half"  of  the  boundary  information  is  known,  j  However,  when  applied  to 
the  boundary,  it  provides  an  alternative  analytical  statement  of  the 
boundary- in i t i al  value  problem  (BIVP).  Consider  applying  Equation  38  to 
an  arbitrary  point  on  T: 

Uk^r,t)  =  /  tGik*ti^-,t-)  dA(^  *  j  I  F  ik*'1  i  1  (  — ' t)  dA(?o 
r  r 

+  j  ( Gik*^i 1 ’ t)  dV(x)  +  T^(£r ,t;vc,Uo)  (39) 

Q 

Note  that  the  second  boundary  integration  involves  an  improper  in¬ 
tegral  to  be  interpreted  as  a  Cauchy  principal  value  intergral.  In  general 
this  can  be  written  as: 


j  [F.k*u.](x,t)  dA(x)  =  FS[kNG  u.(£r,t)  +  f  (F.k*u.l(x,t)  dA(x) 


(40) 


s  i  Nr 

where  F^k  is  the  singularity  contribution  and  the  last  integral  must 
be  intepreted  as  a  Cauchy  principal  value  integral.  For  smooth  boundaries 
F^k^G  =  -1/2  6,k  (see  e.g.,  Cole,  Kosloff,  and  Minster,  1978,  Appendix 
A)  and  Equation  39  becomes: 
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I  Uk(^r>t) 


dA(  x) 


*  /  {lGiknil(x,t)  -  |F|k*u,l(x.t) 

r 

+  f  lG.k*+.Hx,t)  dV(x)  +  I^(£r,t;vo,Uo)  (41) 

n 

This  gives  us  a  boundary  integral  equation  (BIE)  statement  of  the 
boundary- in  it ial  value  problem.  Relations  similar  to  Love's  integral 
identity  can  be  obtained  for  strain,  stress,  and  traction  by  applying 
Equations  2  and  3,  and  Cauchy's  formula,  to  Equation  38;  in  a  numerical 
setting  these  suggested  relations  are  often  not  used.  The  numerical 
approximation  of  the  above  integral  equations  will  be  considered  in 
subsequent  sections.  First,  we  consider  the  integral  equations  which 
the  indirect  boundary  element  method  (1BEM)  approximates  in  the  next 
section. 

Integral  Equations  for  the  Indirect  Boundary  Element  Method.  The 

integral  equations  which  the  IBEM  numerically  approximates  can  be  deduced 
from  Equation  38.  Banerjee  and  Butterfield  (1981)  present  the  analogous 
development  for  steady-state  and  transient  potential  flow  using  an  idea 
originally  due  to  Lamb  (1932).  Eringen  and  Suhubi  (1975,  Section  5.14) 
derive  the  integral  equations  for  elastodynamics  by  the  same  argument. 

An  overview  of  the  derivation  is  presented  below. 

The  indirect  formulation  can  physically  be  visualized  as  embedding 
the  domain  of  the  problem  in  an  infinite  space  of  the  same  media.  The 
traction  distribution,  along  the  surface  corresponding  to  the  boundary 
of  the  original  problem,  is  then  sought  which  will  satisfy  the  boundary 
conditions . 

The  integral  equations  for  this  method  can  be  derived  from  Equation 
38  by  considering  two  displacement  boundary- ini t ial  value  problems: 

1.  A  displacement  BIVP  with  body  forces  and  initial  conditions 
defined  on  a  domain  Q  and  specified  boundary  conditions  on  the  boundary 

r. 
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2,  A  displacement  BIVP  in  the  compliment  space  12  having  the  same 
boundary  conditions  but  with  nil  body  forces  and  initial  conditions. 

Q 

For  a  finite  domain  problem  with  boundary  T,  12  and  12  correspond  to 
the  domains  for  the  interior  and  exterior  problems,  respectively.  Applying 
Equation  38  to  both  problems  and  equating  the  two  expressions  (knowing 
the  boundary  conditions  are  equal)  gives  the  relation: 

u.(x,t)  =  f  [Glk*Pk](£,t)  dA(£)  +  j  |G.k**k](£,t)  dV(£) 

r  12 

+  I V(x , t ; v  ,u  )  (42a) 

l  -  -o-o 

where  the  unknown  vector  valued  function  l’k(£,T)  is  the  difference  in 
the  boundary  traction  distributions  for  the  two  problems,  given  as: 

Pk($,T)  =  tk(S,T)  -  t£(£,i)  (S,T)  e  r  x  T+  (42b) 

is  often  referred  to  as  the  fictitious  or  artificial  traction 
aistribution,  it's  artificial  in  the  sense  that  it  is  a  consequence  of 
embedding  the  problem  in  the  infinite  domain  and  has  no  meaning  outside 
of  this  context.  To  obtain  the  above  relationship,  symmetry  of  the  Green's 
function  with  respect  to  its  indices  and  spacial  arguments  was  employed. 

Given  the  expression  for  the  displacement  field  in  Equation  42a, 
the  kinematic  equation  (Equation  2),  constitutive  relation  (Equation  3), 
and  Cauchy's  formula  then  give  the  traction  field  for  a  unit  normal  nj(x) 
as : 


t1(x,t) 


iFik*pk1(£’t)  dA^>  +  /  dv(^) 

12 


+  Ik(^>t^o’yo) 


(43) 


where 


I.(x,t;vo,uo) 


P  j  Vok(S)  Fik(x,l;^  +  uoR(£)  Fik(x,t;i)]  dV(£) 
12 
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Similar  expressions  for  the  stress  and  strain  fields  are  inter¬ 
mediate  steps  in  obtaining  Equation  43,  but  often  these  expressions  are 
not  used  in  the  numerical  formulation.  Tf  the  artificial  traction  dis¬ 
tribution  were  known,  the  response  in  the  domain  would  be  completely 
defined  by  integral  equations  (e.g..  Equation  42  gives  the  displace¬ 
ment  in  the  domain).  To  attain  the  artificial  tractions,  we  bring  the 
response  point  to  the  boundary  and  enforce  the  boundary  conditions. 

To  obtain  the  key  integral  equation  which  the  IBEM  is  based  on,  we 
assumed  a  displacement  BIVP.  In  actuality,  we  apply  the  IBEM  to  the 
fully-mixed  BIVP  and,  thus,  must  consider  traction  boundary  conditions 
also.  Note  that  while  Equation  42  is  regular  upon  integration,  Equation 
43  must  be  interpreted  in  a  Cauchy  principal  value  sense.  The  boundary 
integral  equations  are  then  given  as: 


u.(xr,t)  =  f  [G.k*Pk](£,t)  dA(£)  +  j  [G.k**k](i,t)  dV(£) 
r  R 

+  I^(xr,t;v  ,u  )  Xp  £  Ty 

;.(xr,t)  =  ij.ik  Pk(xr,t)  +  /  [Fik*PkKi,t)  dA(i) 

r 

+  f  (fik^k](i,t)  dV(£) 

Q 

+  lT(Xp ,t;vo,uo)  xr  £  rT 


(44a) 


(44b) 


where  "  "  and  have  been  used  to  explicitly  indicate  the  known  and 
unknown  field  variables,  respectively.  In  Equation  44b,  the  singular 
contribution  is  based  upon  a  smooth  boundary  and  the  sign  depends  upon 
the  orientation  of  the  normal  vector. 

This  integral  equation  approach  is  sometimes  referred  to  as  "an 
integral  equation  representation  by  vector  simple-layer  potentials," 
which  reflects  its  potential  theory  origins  (see  e.g.,  Jaswon  and  Symm, 
1977).  Equations  44a  and  44b  are  vector  Fredholm  integral  equations  of 
the  first  and  second  kinds,  respectively. 
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Among  the  differences  between  the  integral  equations  on  which  the 
DBF.M  and  IBF.M  are  based,  note  that  for  the  TBFM  integration  of  boundary 
displacements  has  been  eliminated  and  thus  only  a  single  boundary  inte¬ 
gration  in  terms  of  the  artificial  tractions  remains.  Because  of  the 
symmetry  in  the  Green's  tensor,  this  provides  a  field-source  interpre¬ 
tation  to  the  integral  equations  --  a  physically  meaningful  interpreta- 
t  ion . 

In  the  following  sections  we  consider  the  numerical  approximation 
of  the  integral  equations  developed  above  and  how  these  numerical  approxi¬ 
mations  might  be  combined  with  those  of  the  FF,M. 


NUMERICAL  SOLUTION 

The  above  integral  equation  formulations  provide  a  potentially 
effective  approach  to  modeling  infinite  domains.  In  the  first  sub¬ 
sections  below,  we  will  consider  the  numerical  approximations  made  in 
the  integral  equations  to  obtain  the  direct  and  indirect  BEMs .  The 
third  subsection  addresses  how  the  efficiency  of  both  BEMs  can  be 
"optimized"  by  exploiting  properties  of  the  free-space  Green's  function 
and  modern  computer  hardware.  The  final  subsection  gives  an  overview  of 
coupling  BEMs  with  the  FEM  by  treating  the  BFM  as  a  nonlinear  boundary 
condition  on  the  FEM  subregion  of  the  problem. 

Approximation  of  Integral  Equations  --  Boundary  Element  Methods 

The  following  approximations  are  common  to  both  BEMs:  (1)  the 
integrations  are  performed  in  a  piecewise  manner,  and  (2)  the  boundary 
integral  equations  are  approximately  satisfied  in  a  boundary  weighted 
residual  sense  (usually  by  collocation).  Though  by  tradition  we  appear 
to  be  stuck  with  the  name  "boundary  element  methods,"  it  is  somewhat  of 
a  misnomer;  it  suggests  that  analagous  to  finite  elements  (on  the  domain) 
we  will  have  boundary  elements  with  known  shape  functions  which  constrain 
the  boundary  displacement.  We  often  use  the  same  locally  supported  family 
of  polynomials  in  the  BF,M  as  the  FF.M  uses  as  shape  functions  (Lac-hat  and 
Watson,  1975);  however,  these  functions  merely  facilitate  the  piecewise 
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integrations  and  do  not  imply  that  the  fields  are  actually  constrained 
to  these  interpolation  functions.  This  implies  that  a  coupling  between 
the  BE  and  FE  methods  will  be  inherently  incompatible.  For  a  numerical 
demonstration  of  this  incompatibility  in  e 1 astostatics  see  Cox  (1988). 

Let's  consider  the  interpolation  of  boundary  and  domain  quantities 
before  addressing  a  specific  BEM.  (The  presentation  given  here  is  intended 
for  individuals  familiar  with  the  discretizations  commonly  made  in  BEHs; 
for  more  detail,  see  Banerjee  and  Butterfield,  1981  or  Brebbia,  Telles, 
and  Wrobel,  1984.)  The  approximations  for  the  displacement  and  traction 
in  terms  of  expansion  functions  in  time  and  space  (on  T)  are  given  as: 


N  N 
NP  TS 


ui(xr,t)  =  \  \  Nn(~T)  ,]Ij'S 


n=l  s=l 


nnp  nts 


Xj,  e  r 


V-pt)  }  }  Nn(-r)  Ts(t)  ti 


n=l  s=l 


(45a) 


(45b) 


Similarly  if  we  let  x^(Xjj.t)  represent  any  of  the  vector-valued 
field  quantities  in  Equations  38  through  44,  it  can  be  approximated  as: 


nmp  nts 


xi(sa-t)  =  l  l 


where 


-ns 


Mm(x0)  T  (t)  xmS 
m  M  s  i 


X^  €  fi 


(45c) 


m=l  s= 1 


th 


N  (x_)  =  n  expansion  function  for  traction  or  displacement 
n  i  ji 

on  r 


th 


M  (x^)  =  m  expansion  function  for  a  field  variable  on  fi 
m  jc 


th 


T  (t)  =  s  expansion  function  for  time 


N^p  =  number  of  boundary  expansion  functions 


N^p  =  number  of  domain  expansion  functions 


N^,g  =  number  of  time  expansion  functions 


=  expansion  coefficients  for  the  corresponding 
quantities 


-ns  -ms 
i  ’  i 


In  the  more  general  case,  each  boundary  and  domain  expansion  could 

be  written  in  terms  of  unique  expansion  functions.  For  implementation 

purposes,  let's  be  more  specific  instead  of  more  general.  The  time  do- 

A 

main  of  the  analysis  [0,T  ]  is  subdivided  into  N  time  intervals  where 

th  ^ 

the  time  at  the  end  of  the  j  interval  is  denoted  as  t  ^ .  In  general, 
it  is  not  necessary  that  these  time  intervals  (or  steps)  be  of  uniform 
duration;  however,  as  we  will  see  in  a  later  section,  uniform  steps  allow 
a  more  efficient  time-stepping  procedure.  The  boundary  is  discretized 
into  N_.r  surfaces  which  (despite  the  misnomer)  will  be  referred  to  as 
elements;  the  union  of  these  surfaces  spans  the  boundary.  Associated 
with  the  elements  are  nodal  points  occurring  at  the  extremeties  or 
within  the  elements.  Similarly,  the  domain  is  discretized  into  N 
volume  cells;  the  union  of  which  spans  the  portions  of  the  domain  where 
nonzero  body  forces  or  initial  conditions  occur.  Associated  with  the 
cells  are  N^p  cell  points  occurring  on  the  verticies,  edges,  and  within 
the  cells.  Geometrically,  the  boundary  elements  and  volume  cells  are 
similar  to  FE  shell  and  solid  (e.g.,  brick)  elements.  So  that  the  ex¬ 
pansion  coefficients  will  correspond  to  the  nodal  values  of  the  corres¬ 
ponding  boundary  values,  the  following  is  required: 


N  (x.)  =  «  . 

n  ~J  ij 


T  (t.)  =  6  . 

s  J  sj 


x .  e  r 

-j 


(46a) 

(46b) 


.  th 


where  Xj  =  position  vector  for  the  j  nodal  point.  Additionally,  for 
the  expansion  coefficients  of  the  domain  to  correspond  to  the  values  at 
the  cell  points,  we  require: 


M  (x.)  =  6  . 

m  -J  mj 


X  .  E  ft 

“J 


(46c) 


.  th 


where  x.  =  position  vector  for  the  j  "  cell  point. 

—ns 

So,  as  an  example,  u^  of  Equation  45a  is  the  displacement  vector 
at  node  n  and  time  step  s.  As  suggested  at  the  beginning  of  this  section, 
the  interpolation  functions  in  space  often  correpond  to  the  polynomial 
shape  functions  of  the  FEM.  Before  proceeding,  let's  discuss  a  few  im¬ 
portant  exceptions  which  clarify  the  terminology: 
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1.  An  exception  to  this  occurs  with  so-called  "discontinuous  elements." 
The  simplest  example  is  the  constant  element  where  the  integration  over 

the  element  treats  the  boundary  value  as  a  constant.  In  this  case  the 
nodal  point  resides  at  the  center  of  the  element.  Note  that  the  use  of 
the  term  "node"  here  corresponds  to  a  point  where  boundary  values  are 
specified  and,  thus,  for  the  constant  element  is  not  sufficient  to  also 
define  the  geometry.  With  "discont inuous  elements"  it  is  easy  to  define 
loading  discontinuities  if  they  occur  along  element  boundaries;  adjoining 
elements  simply  have  distinct  node  points  approaching  the  element  boundary. 

In  practice  we  let  these  points  coincide  and  simply  perform  the  necessary 
bookkeeping  to  associate  the  corresponding  boundary  values  with  the  correct 
element.  "Discontinuous  elements"  will  be  further  discussed  in  the  next 
section  with  regard  to  the  collocation  technique. 

2.  Another  exception  occurs  for  infinite  domains  where  an  interface 
between  BE  regions  extends  to  infinity,  special  interpolation  functions 
are  required  to  formulate  so-called  "infinite  boundary  elements." 


All  of  these  interpolation  functions  have  the  mathematical  attribute 
referred  to  as  "local  support;"  though  they  are  defined  over  the  whole 
domain/boundary,  they  are  nonzero  only  in  adjoining  cells/elements. 

Local  support  will  not  provide  a  sparse  system  of  equations  as  with  the 
FEM,  however  it  does  reduce  assembly  effort  in  forming  the  systems  of 
algebraic  equations  that  approximate  the  integral  equations. 

The  expansion  functions  in  time  can  also  be  thought  of  as  having 
local  support.  These  functions  are  defined  such  that: 
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That  is,  they  are  polynomial  functions  on  the  time  step  preceeding 
t  .  Different  choices  of  the  polynomial  coefficients  give  various  inte¬ 
gration  rules  in  time. 

These  interpolations  are  the  first  of  two  (principal)  approximations 
to  be  made.  The  second  approximation  addresses  how  the  unknown  boundary 
values  are  obtained;  how  do  we  obtain  the  "best"  expansion  coefficients? 


Direct  Boundary  Element  Method.  This  section  addresses  the  numerical 
approximation  of  integral  equations  38  and  41  --  the  direct  boundary 
element  method.  As  a  stand-alone  analysis  tool  or  in  a  coupled  solution, 
Equation  41  must  be  approximated  on  the  boundary  to  attain  the  unknown 
boundary  values.  In  this  section  we  consider  the  application  of  the 
DBEM  as  a  stand-alone  analysis  tool;  many  of  the  equations  are  common  to 
the  coupled  solution  approach.  The  response  within  the  domain  at  any 
point  can  then  be  calculated  by  integral  equations  like  Equation  38. 

Using  the  interpolations  of  the  previous  section,  the  boundary  integral 
equation  (Equation  41)  is  approximated  an: 


NP  j 

I  -  1  }  r  -r) 

n=l  s=l 


where 


N, 


MP  j 


MP 


+  l  l  *?  +  l 


m=l  s=l 


,  ~m  ,  ,  m 
+  "oi 


m=l 


A  te(tj_1,tj]  (47a) 


g^Ci.t)  =  j  J  Glk(x,t;£,T)  N}(x)  Ts(x)  dt  dA(x) 


r  0 


(47b) 
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t 

fik(£’t}  =  /  f  Fik0i.t-.£.'0  ^(x)  ts(t)  dT  dA(x) 

r  o 

t 

'g“(£,t)  =  f  f  G.k(x,t;i,T)  Mm(x)  Ts(x)  di  dV(x) 

n  o 


~m 

8ik 


a,t) 


f  G.k(x,t;£)  Mm(x)  dVCx) 


(47c) 


(47d) 


(47e) 


iik(i,t)  =  I  G.k(x,t;£)  Mm(x)  dV(x)  (47f) 

Q 

The  local  support  of  each  interpolation  function  tacitly  implies 
reduced  limits  of  integration  in  both  time  and  space. 

The  nodal  boundary  values  at  the  time  t^  are  approximated  by 
satisfying  the  integral  equation  in  a  collocation  sense  with  respect  to 
both  time  and  space.  Other  boundary  weighted  residual  techniques  could 
be  used  but  collocation  is  the  most  prevalent.  In  time,  collocation 
occurs  at  the  ends  of  the  previously  defined  time  steps;  in  space,  the 
nodes  are  the  most  common  collocation  points.  Some  researchers  (see 
e.g.,  Patterson  and  Sheikh,  1981  or  Brebbia,  Telles  and  Wrobel,  1984) 
only  collocate  inside  the  elements  to:  (1)  eliminate  special  singular 
contribution  calculations  when  geometric  discontinuities  coincide  with 
nodes,  and  (2)  to  simplify  the  element  assembly  procedures.  Internally 
collocated  elements  are  often  referred  to  as  "discontinuous"  or  "non- 
conforming"  elements  --  another  misnomer  based  on  the  shape  function 
fallacy.  The  following  development  does  not  exclude  nodal  collocation; 
however,  the  singular  contribution  used  to  obtain  Equation  41  would 
change  for  a  boundary  point  which  does  not  have  a  unique  tangent  plane. 

Applying  Equation  47  to  Np  collocation  points  at  the  end  of  the 
j  time  step  gives: 
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where 
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Note  that  the  vector  gives  the  effect  of  body  forces  and  initial 
conditions  on  the  displacement  at  the  collocation  point.  By  combining 
the  singular  contribution  term  on  the  left-hand  side  with  the  appropriate 

■/r 

DF  matrices,  one  obtains: 


NP  j 

s=  i  Ik 

n=l  s=l 


S  tS  -  DFJS  uS  +  RJ  c  =  l,2,...,N_D 
cn  — n  cn  “n|  c  CP 


th 


(49) 


If  the  collocation  point  coincides  with  the  p  node,  the  relationship 

■Ar 

between  DF  and  DF  is  given  as: 


DFjS  =  *DFjS  +  |  I  6  .6 
cn  cn  2  sj  np 
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3 

where  I  is  the  identity  matrix  for  R  .  Combining  the  equations  associated 
with  each  collocation  point  into  a  single  system  of  equations  gives: 


Assuming  N^p  =  N^p,  we  have  GS,  Fs  t  R  ^  and  tS ,  uS ,  i  R  ^ 
for  all  s  where  =  3N^,p  =  3N^p.  As  an  alternative  one  could  over¬ 
collocate  the  integral  equations  (Ncp>NNp,  see  e.g.,  Hutchinson,  1985) 
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and  then  solve  the  resulting  overdetermined  system  in  a  least  squares 
sense;  however,  its  unlikely  the  overcollocation  would  be  more  effective 
than  solving  the  integral  equations  in  a  Galerkin  sense. 

As  explicitly  denoted  by  the  -'s,  the  only  unknowns  are  half  the 
boundary  values  at  the  current  time  step;  actually  this  is  an  induction 
assumption  which  we  know  is  valid  for  j=t.  In  the  following  equation, 
the  induction  is  validated  and  a  time  stepping  procedure  is  given. 
Partitioning  the  boundary  value  vectors  for  the  time  step  as: 


0  = 


j-1 

l  ( 

s— 1 


GS  tS  -  FS  uS 


)  + 

G^  G^ 

F-!  F;j 

uj 

1 

.  1  2_ 

tj 

i  2 

ij. 

+  RJ 


(51) 


and  isolating  the  unknowns  on  the  left-hand  side  gives  the  following  set 
of  equations  for  the  unknown  boundary  values: 


GJ  -  FJ 

tj 

FJ  -  GJ 

V 

j-1 

]  ( 

1  2. 

uj 

1  2 

tj 

L  \ 

S=1 

GS  tS  -  FS 


uSJ  -  Rj 


(52) 


The  left-hand  side  now  consists  of  the  product  of  a  fully-populated 
coefficient  matrix  with  the  vector  of  unknown  boundary  values  at  time 
t j .  The  right-hand  side  of  the  system  of  equations  is  determined  by 
four  physically  meaningfull  sets  of  values:  (1)  the  known  boundary 
values  at  time  t  ,  (2)  all  the  boundary  values  (displacements  and  trac¬ 
tions)  at  previous  time  steps,  (3)  the  time  variation  of  the  body  loads 
to  time  tj,  and  (4)  the  initial  conditions. 

The  unknown  boundary  values  at  time  L  ^  are  obtained  by  solving  the 

above  system  of  linear  equations.  With  an  approximation  of  the  boundary 

values,  the  response  at  any  point  in  the  domain  over  the  period  of  the 
th 

j  time  step  can  be  determined.  Equation  38  gives  the  necessary  integral 
equation  for  the  displacement  field.  The  numerical  approximation  of 
Equation  38  parallels  the  numerical  approximation  of  Equation  41  except, 
instead  of  considering  a  boundary  collocation  point,  we  consider  an  in¬ 
ternal  (domain)  response  point.  The  following  revisions  of  Equation  47 
yield  the  numerical  approximation  of  Equat  ion  38:  (1)  eliminate  .the  1/2 

factor  on  the  left-hand  side,  and  (2)  interpret  Equation  47c  as  a  regular 
integration  since  £  e  ft. 
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As  noted  in  a  previous  section,  similar  relations  for  the  stress 
and  strain  fields  can  be  easily  derived.  If  one  wanted  both  the  stress 
and  displacement  fields  in  a  region,  it  would  probably  be  more  efficient 
to:  (1)  calculate  the  displacement  at  a  regular  grid  of  points  in  the 
domain,  and  then  (2)  approximate  the  stresses  by  finite  difference  or 
finite  element  approximations.  In  the  later  case  the  derivatives  operate 
on  the  interpolation  functions  which  span  regions  between  points  of  known 
displacement . 

Equation  52  reflects  the  numerical  burden  that  is  characteristic  of 
formulations  which  use  integral  equations  in  time  --  convolution.  The 
following  section  will  discuss  some  ideas  which  significantly  lighten 
this  burden.  These  ideas  on  improving  the  numerical  efficiency  are 
equally  applicable  to  the  indirect  BEM,  so  let's  first  consider  its 
formulation  in  the  following  section. 

Indirect  Boundary  Element  Method.  This  section  addresses  the 
numerical  approximation  of  integral  Equations  42  through  44  --  the  in¬ 
direct  boundary  element  method.  As  a  stand-alone  analysis  tool  or  in  a 
coupled  solution,  Equations  42  and  43  must  be  approximated  on  the  boun¬ 
dary  to  attain  the  artificial  boundary  tractions.  In  this  section  we 
consider  the  application  of  the  IBEM  as  a  stand-alone  analysis  tool; 
many  of  the  equations  are  common  to  the  coupled  solution  approach.  In  a 
coupled  solution  approach,  the  equations  "are  applied  a  second  time"  to 
determine  the  unknown  boundary  values.  The  amount  of  numerical  effort 
is  essentially  equal  to  the  direct  method;  however,  the  indirect  method 
calculates  the  artificial  tractions  as  an  intermediate  step  --  thus  its 
name.  With  the  artificial  tractions,  the  response  in  the  domain  at  any 
point  can  be  calculated  by  integral  equations  like  Equation  42a.  Using 
the  interpolations  of  Equation  45,  boundary  integral  Equation  ^4  is 
written  as: 
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r  0 
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a  o 

gjk(*,t)  =  j  G.kCx,t;0  Mm(4)  dV($) 
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(53a) 

(53b) 

(53c) 

(53d) 

(53e) 


Since  the  free  space  Green's  function  is  symmetric  with  respect  to 
its  spacial  arguments  and  indices,  the  above  integrations  are  equivalent 
to  those  of  Equation  47.  In  a  similar  manner.  Equation  44b  is  approximated 
as : 
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where 
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As  previously  mentioned,  the  local  support  of  the  interpolation 
functions  tacitly  implies  reduced  limits  of  integration  in  both  time  and 
space. 

Again  we  will  limit  the  discussion  to  the  collocation  satisfaction 

of  the  boundary  integral  equations.  With  the  1BEM  we  apply  Equation  53a 

to  collocation  points  on  F^  and  apply  Equation  54a  to  collocation  points 

on  iy  As  previously  mentioned,  we  are  actually  addressing  the  fully- 

mixed  BIVP;  thus,  the  notation  T  and  F  is  symbolic.  Applying  Equations 

^  til 

53a  and  54a  to  N  collocation  points  at  the  end  of  the  j  time  step 

u  I 

gives : 
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~  number  of  displacement  boundary  conditions  at  collocation 
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Note  that  the  vectors  RG^  and  RFjJ  give  the  effects  of  the  body  forces 
and  initial  conditions  on  the  displacement  and  traction  boundary  values, 
respectively.  By  combining  the  singular  contribution  term  with  the  appro- 

it 

priate  DF  matrices  one  obtains: 
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If  the  collocation  point  coincides  with  the  p  node,  the  relation- 
* 

ship  DF  and  DF  is  given  as: 


DFjS  =  *DFjs  ±  |  I  6.  6 

cn  cn  2  js  np 

Let  N,,__  denote  the  total  number  of  equations  associated  with  dis- 
UEQ 

placement  boundary  conditions  and  denote  the  total  number  of  equations 

associated  with  traction  boundary  conditions;  that  is,  if  we  indexed 
and  N^,  for  each  collocation  point  we  would  have: 


-  I  n5  ■ 


where  =  NL,,  =  3N..n  =  3N_,D.  Combining  the  equations  associated 

Uhy  Ihy  hy  Nr  Lr 

with  each  collocation  point  into  a  single  system  of  equations  gives: 
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and  G'  and  F  are  assembled  as  shown  for  Equation  50  except  that  the 
order  of  the  submatrices  varies  as  indicated  in  Equations  55  and  56; 
thus,  some  of  the  submatrices  do  not  even  exist.  For  example,  if  col¬ 
location  point  number  3  corresponds  to  a  point  with  traction  boundary 

1  s 

conditions  only  (i.e.,  NT=3  and  N^O)  then:  (1)  DG^n  ^oes  not  exist  for 

is  3x3 

all  j,s,n;  and  (2)  DF^n  E  R  for  all  j,s,n.  So  we  have 

N-tp-vXN  N_t,_.xNLr. 

Gs  t  R  UEQ  EQ  and  F*  t  R  ^  E<1 

such  that  the  combination  of  these  partition  matrices  gives  square 
matrices.  Isolating  the  unknowns  on  the  left-hand  side  gives  the  fol¬ 
lowing  set  of  equations  for  the  artificial  tractions  at  the  the  end  of 
the  jt^1  time  step, 
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As  previously  discussed,  the  IBEM  solves  for  the  artificial  traction 
values  instead  of  the  unknown  boundary  values.  With  the  artificial  trac¬ 
tions  determined,  unknown  responses  on  the  boundary  or  in  the  domain, 
over  the  interval  of  the  j*'*1  time  step,  can  be  calculated.  As  an  example, 
to  calculate  internal  displacements  an  equation  like  Equation  53a  can  be 
written,  however  x  is  not  necessarily  on  r  ;  this  also  could  easily  be 
expressed  in  matrix  notation. 
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Efficient  Implementation.  To  use  the  term  "efficient"  in  addressing 
the  numerical  implementation  of  the  time  domain  boundary  element  method 
is  a  bit  tongue-in-cheek;  Equations  52  and  59  reflect  the  numerical  burden 
that  is  characteristic  of  formulations  which  use  integral  equations  in 
time  --  convolution.  In  this  section,  we  will  see  that  the  "bad  news" 
of  the  previous  two  sections  can  be  softened  considerably  --  but  not 
forgotten.  We  will  first  consider  a  few  properties  of  the  Green's  func¬ 
tion  which  can  be  exploited  and  then  consider  how  modern  trends  in  com¬ 
puter  architecture,  parallel  processing,  lend  themselves  to  convolution. 

Cole,  Kosloff,  and  Minster  (1978)  appear  to  have  been  the  first  to 
present  the  time  domain  DBEM  for  elastodynamics .  The  development  in 
their  initial  work,  for  the  sake  of  a  simplified  explanation,  was  limited 
to  the  two-dimensional  case  of  antiplane  strain  (a  scalar  BIVP) .  Their 
work  contains  an  excellent  discussion  on  how  the  Green's  function  proper¬ 
ties  can  be  exploited  to  improve  numerical  efficiency  and  how  these  prop¬ 
erties  motivate  their  selection  of  interpolation  functions.  We  have 
previously  selected  a  class  of  interpolation  functions  and  are  concen¬ 
trating  on  numerical  efficiency  in  this  section. 

Initially,  it  appears  that  Equations  52  and  59  require  the  calcula- 
2  2 

tion  of  0(N  N  )  coefficients;  or  following  Cole,  Kosloff,  and  Minster's 

A  2  2  *  g  *  g 

notation  0(N'N„_)  discrete  kernals,  DG^S  and  DF^S.  However,  the  Green's 
TS  CP  cn  cn 

function  has  the  property  of  time  translation  which  can  be  expressed  as: 


Gj.  j  (x, t+At ;£,  t+At)  =  G„(x,t;£,T) 

FjjCx.t+Atj^.T+At)  =  F  (x,t;£,t) 


(60) 


If  the  temporal  interpolation  functions  have  the  same  property, 
that  is: 


Wt+kAt)  =  Tn(t) 

then  one  has  for  the  discrete  kernals  that. 
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For  this  property  to  be  fully  exploited.  Equation  61  requires  the 

analysis  to  be  limited  to  uniform  time  steps.  With  uniform  time  steps 

only  the  full  matrices  relating  the  j t*'  time  step  to  the  lSt  time  step 

must  be  calculated.  This  includes  matrices  resulting  from  the  numerical 

approximation  of  both  boundary  and  domain  integrations.  Exploiting  this 

property  then  reduces  the  number  of  discrete  kernal  calculations  to 
2 

0(N  N  p)  assuming  we  are  able  to  save  matrices  calculated  at  previous 

^  ^  t h 

time  steps.  For  the  DBEM  (see  Equation  52),  this  means  at  the  j  time 

s  s  2. 

step  we  must  have  saved  G  and  F  for  s=l,2, .  .  .  ,  j-1  or  2Npr.(j-l)  floating 

EQ2 

point  numbers.  Similarly,  for  the  IBEM  (see  Equation  59)  N  (j-1)  float- 

LQ 

ing  point  numbers  must  be  saved  for  a  BE  analysis.  For  a  coupled  solution 

2 

approach,  the  IBEM  must  also  save  2NPr.(j-l)  floating  point  numbers  since 

EQ 

both  the  interface  displacements  and  tractions  are  unknown.  This  will 
be  more  apparent  in  the  following  section. 

In  a  coupled  solution  approach  adaptive  adjustment  of  the  time  steps 
in  an  analysis  will  greatly  increase  the  cost  associated  with  the  BE 
subdomain.  If  the  time  step  s '?e  is  only  "occasionally"  changed  the 
piecewise  uniformity  can  be  exploited  but  to  a  lesser  advantage.  The 
determination  of  the  "optimum"  amount  of  adaptive  time  step  adjustment 
is  a  subject  for  numerical  parameter  studies. 

The  time  translation  property  also  indicates  the  possibility  of 
numerical  instability.  Cole,  Kosloff,  and  Minster  (1978)  consider  the 
iterative  process  as  being  similar  to  a  finite  difference  method  on  the 
boundary  where  the  difference  molecule  expands  backward  in  time  with 
each  step.  They  show  that  when  a  linear  interpolation  in  time  is  used 
for  the  tractions,  the  process  is  marginally  stable  in  theory  and  they 
also  indicate  it  has  been  found  to  be  unstable  in  practice. 

Two  more  properties  of  the  Green's  function  allow  us  to  reduce  the 
number  of  unique  coefficient  calculations  and  associated  storage.  The 
first  property  is  causality,  which  simply  says  "there  is  no  response  at 
a  given  point  in  the  domain  due  to  an  impulse  load  until  the  dilational 
wave  has  had  time  to  travel  to  that  point."  It  can  also  be  seen  (Equa¬ 
tion  22)  that  the  response  of  the  point  is  again  quiescent  after  the 
shear  wave  passes.  In  general,  the  nonzero  response  in  the  domain  due 
to  a  concentrated  impulse  load  is  bounded  by  two  spheres  traveling  at 
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the  speeds  of  the  dialational  and  distortional  waves,  respectively.  In 
ft  x  T+  the  region  affected  by  the  impulse  load  is  bounded  by  two  hyper¬ 
cones  with  axes  along  the  time  axis  (see  e.g.,  Cole,  Kosloff,  and  Minster, 
1978).  The  causality  and  post  shear  wave  quiescence  can  be  expressed 
as : 

GjL j Cx,  t;^*)  =  0  (t-T)Cp  <  r  (63a) 

and 

G. .(x,t;£,x)  =  0  (t-r)C_  >  r  (63b) 

rj  s 

respectively. 

We  will  inherently  violate  both  of  these  with  the  discrete  kernals. 

As  indicated  by  Cole  et  al.  (1978),  the  use  of  interpolators  with  separable 
space  and  time  dependence  is  noncausal.  Consider  the  spacial  linear 
interpolation  functions  of  the  FEM.  When  the  wave  front  passes  a  given 
node  the  interpolation  function  instantaneously  becomes  nonzero  in  a 
localized  region.  Since  this  region  can  extend  beyond  the  wave  front 
there  is  a  noncausal  behavior  in  the  numerical  approximation;  that  is, 
our  numerical  approximation  results  in  responses  at  points  prior  to  the 
dilational  waves  arrival.  Cole  et  al.  (1978)  indicate  that  they  do  not 
expect  the  errors  due  to  this  effect  to  be  large  if  very  localized  inter¬ 
polators  are  used  unless  wavelengths  comparable  to  the  node  separation 
are  encountered.  So  the  motivation  for  using  higher  order  elements  to 
reduce  the  number  of  degrees  of  freedom  is  opposed  by  the  motivation  to 
reflect  the  causality  property  in  the  numerical  model. 

The  causality  property  and  post  shear  wave  quiescence  of  the  Green's 
function  cause  many  of  the  discrete  kernals  to  be  zero  if  we  impose  re¬ 
strictions  on  the  spacial  and  temporal  interpolation  functions.  Cole  et 
al.  (1978)  give  general  requirements  since  they  are  motivating  their 
selection  of  interpolation  functions.  The  restriction  is  simply  the 
local  support  of  the  interpolation  functions  in  time  and  space;  for  our 
interpolation  functions  this  is  expressed  as: 

N  (x)  =  0  |x-x  I  >  L  (64a) 

n  —  1 - n1 

T  (t)  =0  t  <  t  ,  or  t  >  t  (64b) 

s  s  - 1  s 
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where  L  is  the  maximum  distance  from  x^  to  any  point  on  an  element  ad¬ 
joining  the  nodal  point.  For  interpolation  functions  associated  with  an 


extreme  node, 

L  is 

bounded  by  the 

maximum  of  the 

two  adjoining  element 

lengths . 

With  the 

local  support 

the  discrete  kernals  reflect  the 

causality 

property 

as : 
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This  gives  a  lower  limit  on  which  discrete,  kernals  are  zero;  for 
some  combinations  of  interpolation  function  and  collocation  points  the  L 
are  not  needed  in  the  inequalities  of  Equation  65.  These  relations 
(Equation  65)  can  be  thought  of  as  the  discrete  analogs  to  the  causality 
property  and  post  shear  wave  quiescence  of  the  Green's  function. 

A  numerical  implementation  can  then  use  inequalities  similar  to  the 
above  to  eliminate  many  discrete  kernal  integrations  and  their  subsequent 
storage.  For  an  out  of  core  routine  (i.e.,  where  the  coefficient  matrices 
are  written  to  disk)  one  might  save  a  Boolean  variable  prior  to  saving 
each  discrete  kernal  where  the  Boolean  variable  has  the  value  of  "true" 
only  if  the  discrete  kernal  values  are  nonzero.  This  would  then  only 
require  a  1-bit  read  for  zero  discrete  kernals  as  apposed  to  a  576-bit 
read  (assuming  a  64-bit  real  word). 

Before  we  leave  the  subject  of  causal  ity,  we  should  note  that  Cole 
et  al.  (1978)  gave  a  criteria  for  the  time  step  as  At  <  L/(2Cp)  where  L 
was  the  minimum  element  length.  They  discuss  the  motivation  for  this 
criteria  in  terms  of  the  backward  causality  cone  (see  the  reference  for 
details).  Other  researchers  (e.g.,  Manolis,  1983  and  Karabalis  and  Beskos, 
1984)  suggest  similar  criteria  for  three-dimensional  problems. 

The  last  properties  of  the  Green's  function  which  can  be  exploited 
to  reduce  the  numerical  effort  are  the  special  translational  and  rotational 
symmetries.  The  translational  symmetry  can  be  expressed  as: 

Gij(x+z.t;i+z,T)  =  G.j(x,t;£,T)  (66) 
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In  general,  these  symmetries  are  not  present  in  the  discrete  kernals 
since  the  geometric  descretization  lacks  spacial  regularity;  however, 
for  a  coupled  problem  the  geometry  of  the  interface  is  relatively  arbi¬ 
trary  and  thus  it  can  be  defined  to  exploit  these  symmetries.  The  sim¬ 
plest  interface  geometry  for  a  dry  dock  might  consist  of  a  half  cylinder 
with  hemispherical  ends. 

An  implementation  which  exploits  these  symmetries  would  probably  be 
highly  specialized  being  based  on  a  given  interface  surface  geometry. 

The  analysis  cost  would  be  reduced  since  there  would  be  fewer  numerical 
integrations  and  associated  storage;  some  of  the  integrations  would  be 
identical  to  previous  integrations  and  some  could  be  obtained  by  the 
orthogonal  transformation  of  previous  integrations.  The  analysis  cost 
would  be  reduced  but  the  implementation  cost  would  be  increased.  Sig¬ 
nificant  bookkeeping  would  be  necessary  to  relate  integrations  to  previous 
integration  results. 

The  last  item  to  consider  in  this  section  is  parallel  processing. 

The  equations  which  arise  in  the  numerical  approximation  of  integral 
equations  lend  themselves  to  parallel  calculations.  Consider  the  forms 
of  Equations  52  and  59.  At  each  step  the  new  coefficient  matrices  could 
be  calculated  in  a  parallel  manner  assigning  each  processor  to  specific 
elements.  In  addition,  the  matrix  multiplications  associated  with  the 
boundary  values  (or  artificial  tractions)  of  previous  steps  could  be 
assigned  to  different  processors.  The  parallel  calculation  of  the  matrix 
multiplications  necessitates  either:  (1)  a  large  amount  of  dedicated 
memory  for  each  processor,  or  (2)  an  architecture  which  performs  10  in  a 
parallel  manner  also  so  the  system  would  not  be  "bus  limited." 

In  this  discussion,  we  have  tacitly  assumed  that  the  individual 
integrations  (i.e.,  Equations  47,  53,  and  54),  are  efficiently  evaluated. 
We  did  not  address  the  evaluation  of  these  integrations.  For  details  on 
these  calculations  see:  Banerjee  and  Ahmad  (1985);  Banerjee,  Ahmad,  and 
Manolis  (1986);  and  Ahmad  and  Banerjee  (1988). 

Assuming  the  BEM  calculations  are  now  efficient  enough  to  be  of 
practical  use,  we  consider  how  to  combine  the  BEM  with  the  FEM  in  the 
next  section.  For  examples  of  the  application  of  the  time  domain  BEM 
see:  Manolis  (1983  and  1984);  Karabalis  and  Beskos  (1984,  1985,  and 
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1986);  Karabalis,  Spyrakos,  and  Beskos  (1984);  Mansur  and  Brebbia  (1985); 
Banerjee  and  Ahmad  (1985);  Banerjee,  Ahmad,  and  Manolis  (1986);  and  Ahmad 
and  Banerjee  (1988). 

Coupling  the  BEMs  with  the  FEM 

In  this  section  we  consider  coupling  the  direct  and  indirect  BEMs 
to  the  FEM.  The  coupling  approach  is  based  on  an  idea  presented  by  Dr. 
Benjamin  Loret  (1987)  during  a  seminar  at  the  University  of  California, 
Davis.  In  his  presentation  he  gave  the  results  of  an  iterative  coupling 
between  the  FEM  and  an  integral  equation  method  for  a  homogeneous  half- 
plane  in  elastostatics .  The  essence  of  his  work  was  the  treatment  of 
the  BE  subdomain  as  a  nonlinear  boundary  condition  to  the  FE  subdomain. 

We  extend  the  idea  to  time  domain  problems  in  this  section.  Very  brief 
outlines  of  algorithms  for  the  IBEM  and  DBF.M  are  given;  the  FEM  is  not 
presented  in  any  detail. 

As  initially  indicated  in  this  report,  the  motivation  for  coupling 
the  two  methods  is  to  let  each  method  model  the  portions  of  the  domain 
for  which  it  is  best  suited.  In  particular,  we  apply  the  BEM  to  the 
infinite  domain  at  whatever  distance  we  are  willing  to  model  the  media 
as  a  homogeneous,  isotropic,  linear  elastic  material.  The  remainder  of 
the  problem,  structure  and  soil,  is  modeled  by  the  FEM  which  is  well 
suited  to  inhomogeneous,  anisotropic,  inelastic  materials. 

The  problem  is  in  general  nonlinear  since  the  constitutive  law 
governing  the  material  behavior  in  the  FE  region  is  inelastic.  We  then 
add  another  nonlinearity  by  treating  the  BF.  region  simply  as  a  nonlinear 
boundary  condition  to  the  FE  region.  The  main  steps  in  the  coupling  are 
outlined  below: 

1.  Assume  is  fixed  during  the  first  time  step.  (Assumes  the 
body  is  at  rest  at  the  beginning  of  the  analysis.) 

2.  With  the  given  displacement  boundary  conditions  on  Tj,  use  the 
FEM  to  calculate  the  generalized  nodal  forces  along  F^. 


3.  Assuming  the  BE  region  displaces  like  its  interpolation  functions 
along  the  boundary,  calculate  the  traction  distribution  on  from  the 
generalized  nodal  forces. 

4.  With  the  given  traction  distribution  at  the  end  of  the  time 
step,  calculate  the  corresponding  displacements  along 

5.  Repeat  steps  2  through  4  until  the  nodal  forces  and  displace¬ 
ments  converge. 

6.  Repeat  steps  2  through  5  for  each  time  step  using  as  the  initial 
displacement  estimate  the  interface  displacement  obtained  at  the  end  of 
the  previous  time  step. 

The  temporal  aspects  of  the  FEM  are  handled  by  a  finite  difference 
approximation;  the  temporal  aspects  of  the  DEM  are  handled  by  a  numerical 
approximation  of  the  convolution  integral.  The  coupling  of  the  two  methods 
iterates  at  each  time  step  to  approximately  satisfy  equilibrium  and  conti¬ 
nuity  at  their  interface. 

We  consider  each  of  the  BE  formulations  below  in  more  detail.  For 
brevity  assume  the  BE  subdomain's  boundary  consists  of  alone.  This 
is  not  a  limitation  of  the  coupling  method. 


DBEM-FEM  Coupling.  Assumming  traction  boundary  conditions, 

Equation  52  can  be  written  as: 

Fj  uj  =  Gj  t j  +  b  (67a) 

where 


j-1 

b  =  J  (gS  tS  -  FS  uS)  +  Rj  (67b) 

s=l 

and  and  t^  are  the  interface  displacements  and  tractions,  respectively. 

An  overview  of  the  coupling  algorithm  is  presented  below.  The  algorithm 
is  described  in  pseudo-code  with  a  Pascal /Modu la  2  dialect.  Consistent 
with  the  mentioned  languages,  supplemental  comments  which  would  not  represent 
actual  code  are  enclosed  with  (*  *)'s.  The  FF.  calculations  are  simply 
represented  by  a  single  call  to  a  procedure  named  FF._Ca  1  cu  1  at  ions  . 
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FOR  j:=  1  TO  N  DO  (*  for  each  time  step  *) 

I  O 

IF  time  steps  are  uniform  THEN 

FOR  i:=  j  TO  2  BY  -1  DO  (*  not  executed  for  j=l  *) 

G1 : =  G1'1 
F* : =  Fi_1 
END_for_i 
Calculate  G^  and 

IF  j=l  THEN  obtain  LU  factorization  of  F^  END_if 
(*  only  factored  once  *) 

ELSE  (*  only  worst  case  is  shown  --  change  on  time  step  size  each 
step  *) 

FOR  i : =  1  TO  j  DO 

Calculate  G1  and  F1 
END_for_i 

Obtain  L(J  factorization  of  F^ 

(*  factored  each  step  for  worst  case  >v) 

F.ND_if 

Let  u^:  =  *  (*  i.e.,  initially  estimate  the  interface 

displacements  at  t .  by  known  interface  displacements 

.  at  Vi  *> 

Calculate  b  by  Equation  67b  (*  contribution  to  known  vector  not 

dependent  on  step  j  *) 

REPEAT  (*  the  nonlinear  boundary  condition  iteration  *) 

(*  Perform  the  FE  calculations  based  on  the  displacement  BCs 
on  Tj..  Note  that  f~*  denote  the  generalized  nodal  forces 

on  rr  *) 


FF,_Calcu  lat  ions(u^  ,  f' ,  iteration  convergence  criteria,...) 

Calculate  the  traction  distribution  t'  on  for  the  BE 

region  by  assuming  nodal  loads  are  obtained  by  weighting 
the  tractions  by  the  interpolation  functions. 

Calculate  a  new  interface,  displacement  estimate  by  solving 
Equation  67a  for  u' .  Note  that  TAJ  factorization 
of  F''  was  previously  obtained. 

UNTIL  nonlinear  boundary  condition  iteration  has  converged 
F.ND_for_  j  . 
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IBEM- 

■FEM  Coupling.  Assuming 

traction  boundary  conditions, 

Equation 

59  can  be 

written  as : 

FJ  |j 

=  tJ  +  b_ 

(68a) 

_  — t 

where 

j-i 

rf| 

bT  - 

:  -  }  FS  PS  - 

(68b) 

— T 

L 

S  =  1 

CP 

ana  and  are  the  artificial  and  actual  interface  tractions, 
respectively.  With  the  artificial  tractions  obtained  by  the  above 
equation,  the  unknown  displacements  along  the  interface  can  be  obtained 
by  Equation  59  as: 


-J  _  gJ  PJ  +  b 


-U 


where 


j-1 

l 

S=1 


GS  PS  + 


RG; 


RG; 


CP 


(69a) 


(69b) 


These  two  systems  of  equations  can  be  used  to  couple  the  IBEM  to 
the  FEM  analagous  to  the  approach  used  above  for  the  DBEM.  An  overview 
of  the  coupling  algorithm  is  presented  below. 
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FOR  j : =  1  TO  N  DO  (*  for  each  time  step  *) 

I  o 

IF  time  steps  are  uniform  THEN 

FOR  i:=  j  TO  2  BY  -1  DO  (*  not  executed  for  j=l  *) 

Gi:=  Gi_1 
Fi_1 

END_for_i 

Calculate  G^  ana  F^ 

IF  j=l  THEN  obtain  L (I  factorization  of  F^  END_if 
(*  only  factored  once  *) 

ELSE  (*  only  worst  case  is  shown  --  change  on  time  step  size  each 
step  *) 

FOR  i:=  1  TO  j  DO 

Calculate  G1  and  F1 
END_for_i 

Obtain  LU  factorization  of  F^ 

(*  factored  each  step  for  worst  case  *) 

F.ND_it 

Let  u^:=  u^  ^  (*  i.e.,  initially  estimate  the  interface 
displacements  at  t.  by  known  interface 
displacements  at  t ^  ^  *) 

Calculate  b^  and  b^  by  Equations  68b  and  69b,  respectively 

(*  contributions  to  known  vectors  not  dependent  on  step  j  *) 
REPEAT  (*  the  nonlinear  boundary  condition  iteration  *) 

(*  Perform  the  FE  calculations  based  on  the  displacement  BCs 
on  Tj.  Note  that  f'  denote  the  generalized  nodal  forces  on 

rr  ,v)  ..  .. 

FE_Calculations(u^ ,  £ ^ ,  iteration  convergence  criteria,  ...) 
Calculate  the  traction  distribution  f'  on  for  the  BE  region 
by  assuming  nodal  loads  are  obtained  by  weighting  the 
tractions  by  the  interpolation  functions. 

Calculate  a  new  artificial  traction  vector  estimate  by  solving 
Equation  68a  for  P* 

Note  that  LU  factorization  of  F*  was  previously  obtained 
Calculate  a  new  interface  displacement  estimate,  u^ ,  by 
Equation  69a 

UNTIL  nonlinear  boundary  condition  iteration  has  converged 
END  for  i . 


Note  that  tne  algorithms  for  coupling  the  direct  and  indirect  BEMs 
are  very  similar.  For  problems  with  zero  body  forces  and  a  quiescent 
past  the  two  coupled  approaches  require  the  same  amount  of  numerical 
effort;  when  body  forces  or  nonzero  intial  conditions  exist  the  direct 
method  requires  less  computational  effort.  However,  if  many  responses 
within  the  BE  region  are  obtained  in  an  analysis  the  1BEM  could  require 
less  computational  effort  since  it  integrates  the  effect  of  a  single 
time  varying  vector  (the  artificial  tractions)  instead  of  two  vectors 
(the  boundary  displacements  and  tractions)  like  the  DBEM. 

Note  that  the  above  algorithms  do  not  explicitly  show: 

•  exploitation  of  causality  property  and  post  shear  wave 
quiescence 

•  details  of  how  to  exploit  the  time  translation  property  when  the 
time  step  size  changes  only  intermittently 

•  details  of  FE  calculations 

•  internal  response  calculations  for  the  BE  region 

•  parallel  calculation 

•  convergence  criterias 

As  previously  mentioned,  the  displacements  along  the  boundary  of 
the  BE  region  do  not  agree  with  the  interpolat  ion  functions  except  at 
the  collocation  points.  Thus,  the  use  of  the  interpolation  functions  as 
shape  functions,  which  they  are  not,  to  calculate  the  traction  distribution 
along  the  interface  from  the  generalized  nodal  forces  is  inherently  in 
error;  it  also  implies  that  we  only  satisfy  continuity  at  the  collocation 
points.  The  inconsistency  between  >he  assumed  and  actual  d i s t r ibut ion 
of  the  displacements  could  possibly  be  accounted  for  with  a  "weighting 
approach"  by  calculating  the  displacement  at  boundary  points  between 
collocation  points.  At  worst  (assuming  we  have  convergence)  the  incon¬ 
sistency  is  a  good  indicator  of  when  mesh  refinement  is  necessary. 


In  the.  above  algorithms  we  explicitly  showed  the  FEM  calculations 
within  the  iterative  loop  for  the  nonlinear  boundary  conditions  (BCs). 
The  FEM  calculations  themselves  contain  iterations  to  satisfy  the  con¬ 
stitutive  relations  (among  other  things).  Whether  the  constitutive 
relations  should  be  iterated  to  convergence  within  each  BC  iteration  is 
a  subject  for  numerical  studies.  We  are  uncertain  of  the  effects  on 
numerical  efficiency  and  convergence.  The  BEM  calculations  in  each  BC 


smaller  subregions);  thus,  in  the  above  algorithms  we  sought  to  minimize 
the  number  of  BC  iterations.  It  might  be  most  effective  to  have  the 
convergence  criteria  of  the  FE  calculations  become  tighter  as  the  BC 
iterations  converge. 


CONCLUSIONS 

The  time  domain  boundary  element  methods  based  upon  the  Stokes 
solution  appear  to  be  the  best  suited  BE  formulations  for  the  coupled 
solution  of  structural/geotechnical  interaction  problems  that  include 
nonlinearities  and  infinite  domains.  The  direct  boundary  element  method 
(DBEM)  has  a  very  elegant  analytical  basis  which  expresses  the  elasto- 
dynamic  boundary-initial  value  problem  in  an  integral  equation  form  using 
the  Stokes  solution  and  the  dynamic  equivalent  to  Betti's  reciprocal 
work  theorem.  The  integral  equations  for  the  indirect  boundary  element 
method  can  be  derived  from  those  for  the  DBEM  by  considering  an  exterior 
and  interior  problem  with  a  common  boundary.  An  artificial  traction 
along  the  common  boundary  is  determined  which  satisfies  the  boundary 
conditions  of  the  actual  problem  while  displacement  compatabi 1 ity  of  the 
two  problems  is  enforced.  The  analytical  basis  for  both  methods  involves 
integral  equations  in  both  time  and  space. 

While  the  integral  equations  are  very  elegant  theoretical  formula¬ 
tions,  their  numerical  approximations  are  computationally  very  demanding. 
For  more  efficient  computations  one  can  exploit  certain  properties  of 
the  Green's  function  and  modern  computer  hardware,  such  as  parallel 
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processing;  this  requires  considerable  algorithm  implementation  effort 
but  appears  to  be  a  necessity  for  anything  but  an  academic  computer 
program. 

Time  domain  BEMs  based  upon  Stokes’s  solution  provide  a  very  rig¬ 
orous  solution  to  the  radiation  boundary  condition  problem  associated 
with  infinite  domains.  The  theoretical  and  numerical  background  for  the 
methods  provides  a  good  basis  for  future  implementation  work.  As  to 
whether  the  coupled  solution  approach  wi1!  provide  a  cost-effective 
solution  for  problems  with  infinite  domains,  speculation  is  a  poor  sub¬ 
stitute  for  numerical  experience.  However,  careful  algorithm  design 
and  computer  hardware  advances  will  improve  the  potential  of  applying 
the  methods  to  large  structural/geotechnical  problems. 


RECOMMENDATIONS 

The  following  outlines  the  recommended  steps  for  the  development  of 
a  coupled  FEM/BEM  computer  code  for  the  dynamic  analysis  of  soil-structure 
interaction  problems: 

1.  Investigate  the  numerical  integrations  necessary  for  the  BEM  in 
elastodynamics .  That  is:  (a)  classify  the  different  integrations  (e.g., 
singular  versus  nonsingular),  (b)  determine  appropriate  numerical  integra¬ 
tion  schemes,  and  (c)  investigate  the  accuracy  of  the  numerical  schemes. 

2.  Develop  a  simple  stand-alone  BEM  computer  program  for  elasto¬ 
dynamics  . 

3.  Extend  the  BEM  program  to  exploit  Green's  function  properties. 

4.  Extend  the  BEM  program  to  exploit  parallel  processing  hardware. 

5.  Extend  the  BEM  program  for  use  in  a  coupled  FEM-BEM  solution. 
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NOMENCLATURE 


A  partial  list  of  symbols  used  in  this  report  follows. 

Mathematical  Symbols 


[$*X](x,t) 


0 

t 

l  $(x,t-T) 

0 


X(x,t)  dt, 


(x,  t)  E  X  T 

(x,t)  e  x  T+ 


where  $  and  X  are  such  that  { $(x,t) ,x(x, t) ]  z  C  (fl  x  T  ). 

This  defines  the  Riemann  Convolution.  For  brevity  it  is  often  denoted 
as  simply  #*x. 


Script  Symbols 

E  Class  of  all  elastodynamic  states. 

Eq  Class  of  all  elastodynamic  states  of  quiescent  past. 

y  An  elastodynamic  state. 


Latin  Symbols 


Cm,n(fi  x  T) 


E3 

E3$ 

Rn 


Defines  the  class  of  all  functions  which  have  continuous 
spacial  and  temporal  derivatives  of  order  up  to  and 
including  m  and  n,  respectively,  on  8  x  T. 

, .  3 

Three-dimensional  Euclidian  space  corresponding  to  R  . 
Three-dimensional  Euclidian  space  omitting  the  point  £■ 
The  linear  space  of  ordered  n-tuples  of  real  numbers. 


The  time  interval  [0,«).  Where  [  and  )  denote  the 
closed  and  open  ends  of  the  interval,  respectively. 


T  The  time  interval  (  -«• ,  0 ]  . 

T**  =  T+  U  T" 
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Greek  Symbols 


6  .  . 
ij 

6(x,z) 

t  ,  . 

ij 

r 

0  .  . 

IJ 

a 


*i 


Kronecker  delta  symbol. 

Dirac  'delta  function.' 

Strain  tensor  in  rectangular  Cartesian  system. 

Complete,  finite  boundary  of  the  problem. 

Stress  tensor  in  rectangular  Cartesian  system. 

Domain  of  the  problem. 

Vector  of  body  forces  per  unit  volume. 

Vector  of  body  forces  per  unit  mass. 

Mass  density  at  a  given  point  in  the  domain. 

Position  vector  to  a  point  on  the  boundary  of  the  problem. 
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Island.  Code  280.  Vallejo.  CA.  Mare  Island.  PWO,  Vallejo.  CA:  Norfolk.  Code  380.  Portsmouth.  V'A; 
Norfolk.  Code  440.  Portsmouth.  VA;  PWO.  Bremerton.  WA 
N  A  VST  A  CO.  Long  Beach.  CA;  CO,  Roosevelt  Roads,  PR:  Code  N4214.  Mayport.  FL;  Engr  Div,  PWD. 
Rodman.  Panama  Canal.  Engrg  Dir.  PWD.  Rota.  Spain.  SCE.  San  Diego.  CA;  WC  43.  Guantanamo  Bay. 
Cuba 

NAVSUPPACT  PWO.  Naples.  Italy 

NAVSCPPFAC  Contract  Admin  Lech  Library.  Diego  Garcia 
NAVSCPSYSCOM  Code  0622.  Washington.  DC 

NAVSWC  Code  E2I1  (Miller).  Dahlgren.  VA;  DEL.  While  Oak  Lab.  Code  WSO.  Silver  Spring.  MD;  W4ICI. 
Dahlgren.  V'A 

NAVWARCOI.  Code  24.  Newport.  RI 

NAVWPNCEN  AROICC.  China  Lake.  CA:  PWO  (Code  266).  China  Lake.  CA 

NAVWPNSIA  Code  042B  (Hunt).  Yorktown.  VA.  Dir.  Maint  Control.  PWD.  Concord.  CA;  PWO.  Charleston 
SC;  PWO.  Seal  Beach.  CA;  PWO.  Yorktown.  VA 
NAVWPNSL'PPCEN  PWO.  Crane.  IN 
NEK  Code  42.  Newport.  RI;  PWO.  Newport,  RI 
NCR  20.  CO 

NF.ESA  Code  1IIE  (McClamc).  Port  Hueneme.  CA 
NMCB  3.  Ops  Offr;  40.  CO;  5.  Ops  Dept 
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NCAA  Joseph  Vadus.  Rockville.  MD 

NRL  Code  2511.  Washington.  DC;  Code  4670  (B  Faraday).  Washington.  DC 
NSC  Code  54.1.  Norfolk.  VA 
NSD  SCE.  Suhic  Bay.  RP 

NUSC  DET  Code  2145  (Varlev).  New  London.  CT;  Code  44  (RS  Munn).  New  London.  CT;  Code  TA151,  New 
London.  CT:  Lib  (Code  4555).  Newport.  Rl 
PACMISRANFAC  HI  Area.  PWO.  Kekaha.  HI 

PHIBCB  1.  CO.  San  Diego,  CA;  1.  P&L.  San  Diego.  CA:  2.  CO.  Norfolk.  VA 
PMTC  Code  1018,  Point  Nlugu.  CA;  One  5041.  Point  Mugu.  CA 

PWC  ACE  Office.  Norfolk,  VA:  Code  1  ).  Great  Lakes.  IL;  Code  10.  Oakland.  CA;  Code  101  (Library). 
Oakland.  CA:  Code  1011,  Pearl  Harbor.  HI;  Code  102.  Oakland.  CA;  Code  125-C.  San  Diego.  CA:  Code 
50.  Norfolk.  VA;  Code  400.  Great  Lakes.  IL;  Code  400.  Oakland.  CA:  Code  4(H).  Pearl  Harbor.  HI:  Code 
400.  San  Diego.  CA;  Code  420,  Great  Lakes.  IL:  Code  420.  Oakland.  CA:  Code  420B  (Waid).  Subic  Bay. 
RP:  Code  421  (Kaya).  Pearl  Harbor.  HI;  Code  421  (Ouin).  San  Diego.  CA:  Code  421  (Reynolds).  San 
Diego.  C'A;  Code  422.  San  Diego.  CA;  Code  425,  San  Diego.  CA:  Code  42'  Norfolk.  VA;  Code  450  (Kyi). 
Pearl  Harbor,  HI;  Code  500.  Great  Lakes.  IL;  Code  500.  Oakland.  CA:  Library  (Code  154).  Pearl  Harbor. 
HI;  Library.  Guam.  Mariana  Islands:  Library.  Norfolk.  VA;  Library.  Pensacola.  FL;  Library.  Yokosuka. 
Japan;  Tech  Library.  Suhic  Bay.  RP 
SPCC  PWO  (Code  08X1.  Mechanicshurg.  PA 
SUBASE  Bangor.  PWO  (Code  8525).  Bremerton.  WA 
SUPSHIP  Tech  Library.  Newport  News.  VA 

USDA  Ext  Serv  (T  Maher),  Washington.  DC:  For  Sve  Reg  8.  (Bowers).  Atlanta.  GA;  For  Svc.  Reg  Bridge 
Engr.  Aloha.  OR;  For  Svc.  Tech  Engrs.  Washington.  DC 
USNA  C'h.  Mech  Engrg  Dept.  Annapolis,  MD;  Ocean  Engrg  Dept  (McCormick).  Annapolis.  MD;  PWO. 
Annapolis,  MD 

CALIFORNIA  STATE  UNIVERSITY  C'.V  Chelapati.  Long  Beach.  CA 
CASE  WESTERN  RESERVE  UNIV  CE  Dept  ( Perdikaris).  Cleveland.  OH 
CATHOLIC  UNIV  of  Am.  CE  Dept  (Kim).  Washington.  DC 
CITY  OF  LIVERMORE  Dackins,  PE.  Livermore.  C'A 
CLARKSON  COLL  OF  TECH  C'E  Dept  (Batson).  Potsdam.  NY 
COLORADO  STATE  UNIVERSITY  CE  Dept  (Criswell).  Ft  Collins.  CO 

CORNELL  UNIVERSITY  Civil  &  Environ  Engrg  (Dr  Kulhawy),  Ithaca.  NY;  Library.  Ithaca.  NY 
DAMES  &  MOORE  Library.  Los  Angeles.  CA 

FLORIDA  ATLANTIC  UNIVERSITY  Ocean  Engrg  Dept  (Su).  Boca  Raton.  FL 
FLORIDA  INST  OF  TECH  CE  Dept  (Kalajian).  Melbourne.  FL 

GEORGIA  INSTITUTE  OF  TECHNOLOGY  CF,  Scol  (Kahn).  Atlanta.  GA:  CE  Scol  (Swanger),  Atlanta,  GA; 
CE  Scol  (Zuruck).  Atlanta.  GA 

INSTITUTE  OF  MARINE  SCIENCES  Library.  Port  Aransas.  IX 
JOHNS  HOPKINS  UNIV  CE  Dept  (Jones).  Baltimore.  MD 

LAWRENCE  LIVERMORE  NATL  LAB  FJ  Tokarz,  Livermore.  CA;  Plant  Engrg  Lib  (L-654),  Livermore,  CA 

LEHIGH  UNIVERSITY  Linderman  Library.  Bethlehem.  PA 

LONG  BEACH  PORT  Engrg  Dir  (Allen),  Long  Beach.  CA 

MICHIGAN  TECH  UNIVERSITY  CE  Dept  (Haas).  Houghton.  MI 

MIT  Engrg  Lib.  Cambridge.  MA;  Lib,  Tech  Reports.  Cambridge.  MA 

NATL  ACADEMY  OF  SCIENCES  NRC.  Naval  Studies  Bd.  Washington.  DC 

NEW  MEXICO  SOLAR  ENERGY  INST  Dr.  Zwibel.  Las  Cruces.  NM 

OREGON  STATE  UNIVERSITY  CE  Dept  (Hicks).  Corvallis.  OR 

PENNSYLVANIA  STATE  UNIVERSITY  Gotolski,  University  Park.  PA;  Rsch  Lab  (Snyder),  State  College. 

PA 

PORTLAND  STATE  UNIVERSITY  Engrg  Dept  (Migliori).  Portland.  OR 

PURDUE  UNIVERSITY  CE  Scol  (Leonards).  W.  Lafayette.  IN;  Engrg  Lib.  W  Lafayette.  IN 

SAN  DIEGO  PORT  Port  Fac,  Proj  Engr,  San  Diego.  CA 

SAN  DIEGO  STATE  UNIV  C'E  Dept  (Krishnamoorthy),  San  Diego.  CA 

SEATTLE  PORT  W  Ritchie.  Seattle.  WA 

SEATTLE  UNIVERSITY  CE  Dept  (Schwacgler),  Seattle,  WA 

SOUTHWEST  RSCH  INST  Energetic  Sys  Dept  (Esparza).  San  Antonio,  TX;  King.  San  Antonio.  TX;  M 
Polcyn.  San  Antonio.  TX;  Marchand,  San  Antonio,  TX 
STATE  UNIVERSITY  OF  NEW  YORK  CE  Dept  (Reinhorn).  Buffalo.  NY;  CE  Dept.  Buffalo.  NY 
TEXAS  A&M  UNIVERSITY  CE  Dept  (Machcmehl),  College  Station.  TX;  CE  Dept  (Niedzwecki),  College 
Station.  TX,  Ocean  Engr  Prop  College  Station,  TX 

UNIVERSITY  OF  CALIFORNIA  C’E  Dept  (Fenves),  Berkeley.  CA;  CE  Dept  (Fourney),  Los  Angeles,  CA; 

CE  Dept  (Gcrwick>,  Berkeley,  CA;  CE  Dept  (Taylor).  Davis,  CA;  CE  Dept  (Williamson).  Berkeley.  CA; 
Naval  Archt  Dept.  Berkeley,  CA 

UNIVERSITY  OF  HARTFORD  CE  Dept  (Kcshawarz).  West  Hartford.  CT 

UNIVERSITY  OF  HAWAII  CF.  Dept  (Chiu),  Honolulu,  HI;  Manoa,  Library.  Honolulu.  HI,  Ocean  Engrg 
Dept  (Ertckin).  Honolulu,  HI 

UNIVERSITY  OF  ILLINOIS  Library.  Urbana.  IL;  Metz  Ref  Rm.  Urbana.  IL 
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UNIVERSITY  OF  MICHIGAN  CE  Dept  (Richart).  Ann  Arbor.  MI 
UNIVERSITY'  OF  NEBRASKA  Polar  Ice  Coring  Office.  Lincoln.  NE 

UNIVERSITY  OF  NEW  MEXICO  HL  Schreyer.  Albuquerque.  NM;  NMERI  (Bean).  Albuquerque.  NM; 

NMERI  (Falk),  Albuquerque.  NM:  NMERI  (Leigh).  Albuquerque.  NM 
UNIVERSITY  OF  PENNSYLVANIA  Dept  of  Arch  (P.  McCleary).  Philadelphia.  PA 
UNIVERSITY  OF  RHODE  ISLAND  CE  Dept  (Kovacs).  Kingston.  RI;  CE  Dept.  Kingston.  RI 
UNIVERSITY’  OF  TEXAS  CE  Dept  (Thompson).  Austin.  TX;  Construction  Industry  Inst.  Austin.  TX;  ECJ 

4.S  (Breen).  Austin.  TX 

UNIVERSITY  OF  WASHINGTON  CE  Dept  (Mattock).  Seattle.  WA 

UNIVERSITY  OF  WISCONSIN  Great  Lakes  Studies  Cen.  Milwaukee.  WI 

WASHINGTON  DHHS.  OFE  PHS  (Ishihara).  Seattle.  WA 

ADVANCED  TECHNOLOGY.  INC  Ops  Cen  Mgr  (Bednar).  Camarillo.  C A 

AMERICAN  CONCRETE  INSTITUTE  Library.  Detroit.  MI 

ARVID  GRANT  &  ASSOC  Olympia.  WA 

ATLANTIC  RICHFIELD  CO  RE  Smith.  Dallas.  TX 

BATTELLE  D  Frink.  Columbus.  OH 

BECHTEL  CIV  IL.  INC  Woolston.  San  Francisco.  CA 

BETHLEHEM  STEEL  CO  Engrg  Dept  (Dismuke).  Bethlehem.  PA 

BRITISH  EMBASSY  Set  &  Tech  Dept  (Wilkins).  Washington.  DC 

BROWN  &  ROOT  Ward,  Houston.  TX 

CANADA  Viateur  De  Champlain.  D.S.A..  Matane.  Quebec 

CHEVRON  OIL  FLD  RSCH  CO  Strickland.  La  Habra.  CA 

CHILDS  ENGRG  CORP  K.M.  Childs.  Jr.  Medfield,  MA 

CLARENCE  R  JONES.  C'ONSULTN.  LTD  Augusta.  GA 

COLLINS  ENGRG.  INC  M  C.arlich.  Chicago.  IL 

CONRAD  ASSOC  Luisoni,  Van  Nuys.  CA 

CONSOER  TOWNSEND  &  ASSOC  Schramm.  Chicago.  IL 

CONSTRUCTION  TECH  LABS.  INC  G.  Corley.  Skokie.  IL 

DAVY  DRAVO  Wright.  Pittsburg.  PA 

DILLINGHAM  CONSTR  CORP  (HD&C),  F  McHale.  Honolulu.  HI 
1  Ong  Yam  Chai.  Singapore 

EARL  &  WRIGHT  CONSULTING  ENGRGS  Jensen.  San  Francisco.  C'A 

EVALUATION  ASSOC.  IN<~  MA  Fedeic.  King  of  Prussia,  PA 

GRUMMAN  AEROSPACE  CORP  Tech  Info  Ctr,  Bethpage.  NY 

HALEY  &  ALDRICH.  INC.  T.C.  Dunn,  Cambridge.  MA 

HARTFORD  STEAM  BOILER  INSP  &  INS  CO  Spinclli.  Hartford.  CT 

HAYNES  &  ASSOC  H.  Haynes.  PE.  Oakland.  CA 

HIRSCH  &  CO  L  Hirsch,  San  Diego.  CA 

HJ  DEGENKOLB  ASSOC  W  Murdough,  San  Francisco.  CA 

HUGHES  AIRCRAFr  CO  Tech  Doc  Cen.  El  Segundo.  CA 

INTL  MARITIME.  INC  D  Walsh.  San  Pedro.  CA 

IRE-ITTD  Input  Proc  Dir  (R.  Danford).  Eagan.  MN 

JOHN  J  MC  MULLEN  ASSOC  Library.  New  York.  NY 

LEO  A  DALY  CO  Honolulu.  HI 

LIN  OFFSHORE  ENGRG  P  Chow.  San  Francisco  CA 

LINDA  HALL  LIBRARY  Doc  Dept.  Kansas  Citv.  MO 

MARATHON  OIL  CO  Gamble.  Houston.  TX 

MARITF.CH  ENGRG  Donoghue,  Austin.  TX 

MC  CLELLAND  ENGRS.  INC  Library.  Houston.  TX 

MOBIL  R&D  CORP  Offshore  Engrg  Lib.  Dallas.  TX 

MT  DAVISSON  CE.  Savoy.  IL 

EDWARD  K  NODA  &  ASSOC  Honolulu.  HI 

NEW  ZEALAND  NZ  Concrete  Rsch  Assoc,  Library.  Porirua 

NUHN  &  ASSOC  A  C.  Nuhn,  Wayzata.  NM 

PACIFIC  MARINE  TECH  (M.  Wagner)  Duvall.  WA 

PILE  BUCK.  INC  Smoot,  Jupiter,  FL 

PMB  ENGRG  Coull.  San  Francisco.  CA 

PORTLAND  CEMENT  ASSOC  AE  Fiorato.  Skokie.  IL 

PRESNELL  ASSOC,  INC  DG  Prcsncll.  Jr,  Louisville,  KY 

SANDIA  LABS  Library.  Livermore.  CA 

SARGENT  &  HERKES,  INC  JP  Pierce.  Jr.  New  Orleans.  LA 

SAUDI  ARABIA  King  Saud  Univ,  Rsch  Cen,  Riyadh 

SEATECH  CORP  I  croni,  Miami,  FL 

SHELL  OIL  CO  E  Doyle.  Houston,  TX 

SIMPSON.  GUMPERTZ  &  HEGER,  INC  E  Hill.  CE.  Arlington.  MA 

TRW  INC  Crawford.  Redondo  Beach.  CA:  Dai.  San  Bernardino,  CA;  Engr  Library.  Cleveland,  OH;  Rodgers. 

Redondo  Beach,  CA 
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TUDOR  ENC.RG  CO  Ellegood.  Phoenix.  AZ 

VSE  Ocean  Engrg  Gp  (Murton).  Alexandria.  VA 

VULCAN  IRON  WORKS.  INC  DC  Warrington.  Chattanooga.  TN 

WESTINGHOUSE  ELECTRIC  CORP  Library,  Pittsburg.  PA 

WISS.  JANNEY.  ELSTNER.  &  ASSOC  DW  Pfeifer,  Northbrook,  IL 

WISWELL.  INC  G.C.  Wiswell.  Southport.  SC 

WOODWARD-CLYDE  CONSULTANTS  West  Reg,  Lib.  Oakland,  CA 

BROWN.  ROBERT  University.  AL 

BULLOCK.  TE  La  Canada.  CA 

CHAO.  JC'  Houston.  TX 

CLARK.  T.  Redding.  CA 

GIORDANO.  A.J.  Sewell.  NJ 

HARDY.  S.P  San  Ramon.  CA 

HAYNES.  B.  No  Stonington.  CT 

HEUZE.  F  Alamo.  CA 

KOSANOWSKY.  S  Pond  Eddy.  NY 

NIEDORODA.  AW  Gainesville.  FL 

PETERSEN.  CAPT  N.W.  Pleasanton.  CA 

QUIRK.  J  Panama  City,  FL 

SPIELVOGEL.  L  Wyncote.  PA 

STEVENS.  TW  Dayton.  OH 

VAN  ALLEN.  B  Kingston.  NY 


68 


INSTRUCTIONS 


The  Naval  Civil  Engineering  Laboratory  has  revised  its  primary  distribution  lists.  The  bottom  of  the 
label  on  the  reverse  side  has  several  numbers  listed.  These  numbers  correspond  to  numbers  assigned  to 
the  list  of  Subject  Categories.  Numbers  on  the  label  corresponding  to  those  on  the  list  indicate  the 
subject  category  and  type  of  documents  you  are  presently  receiving.  If  you  are  satisfied,  throw  this  card 
away  (or  file  it  for  later  reference). 

It  you  want  to  change  what  you  are  presently  receiving: 

•  Delete  -  mark  off  number  on  bottom  of  label. 

•  Add  -  circle  number  on  list. 

•  Remove  my  name  from  all  your  lists  -  check  box  on  list. 

•  Change  my  address  -  line  out  incorrect  line  and  write  in  correction  (DO  NOT  REMOVE  LABEL). 

•  Number  of  copies  should  be  entered  after  the  title  of  the  subject  categories  you  select. 

Fold  on  line  below  and  drop  in  the  mail. 

Note  Numbers  on  label  but  not  listed  on  questionnaire  are  for  NCEL  use  only,  please  ignore  them. 


Naval  Civil  Engineering  Laboratory 
Port  Hueneme.  CA  93043-5003 


Official  Business 

Penalty  for  Private  Use.  $300 


BUSINESS  REPLY  CARD 

FIRST  CLASS  PERMIT  NO  12503  WASH  D  C 
POSTAGE  WILL  BE  PAID  BY  ADDRESSEE 


NO  POSTAGE 
NECESSARY 
IF  MAILED 
IN  THE 

UNITED  STATES 


Commanding  Officer 
Code  L34 

Naval  Civil  Engineering  Laboratory 
Port  Hueneme,  California  93043-5003 


DISTRIBUTION  QUESTIONNAIRE 

The  Naval  Civil  Engineering  Laboratory  is  revising  its  Primary  distribution  lists 


SUBJECT  CATEGORIES 

1  SHORE  FACILITIES 

2  Construction  methods  and  materials  (including  corrosion 

control,  coatings) 

3  Waterfront  structures  (maintenance'deterioratior.  control) 

4  Utilities  (including  power  conditioning) 

5  Explosives  safety 

6  Aviation  Engineering  Test  Facilities 

7  Fire  prevention  and  control 

8  Antenna  technology 

9  Structural  analysis  and  design  (including  numerical  and 

computer  techniques) 

10  Protective  construction  (including  hardened  shelters. 

shock  and  vibration  studies) 

11  Soil  rock  mechanics 

14  Airfields  and  pavements 

15  ADVANCED  BASE  AND  AMPHIBIOUS  FACILITIES 

16  Base  facilities  (including  shelters  power  generation,  water 

supplies) 

17  Expedient  roads/airfields/bridges 

18  Amphibious  operations  (including  oreakwaters  wave  forces) 

19  Over-the-Beach  operations  (including  containerization 

material  transfer  lighterage  and  cranes) 

2u  r-UL  storage  transfer  and  distribution 


TYPES  OF  DOCUMENTS 

85  Techdata  Sheets  86  Technical  Reports  and  Technical  Notes 
83  Table  of  Contents  &  Index  to  TDS 


28  ENERGY/POWER  GENERATION 

29  Thermal  conservation  (thermal  engineering  of  buildings.  HVAC 

systems,  energy  loss  measurement,  power  generation) 

30  Controls  and  electrical  conservation  (electrical  systems. 

energy  monitoring  and  control  systems) 

31  Fuel  flexibility  (liquid  fuels,  coal  utilization  energy 

from  solid  waste) 

32  Alternate  energy  source  (geothermal  power  photovoltaic 

power  systems,  solar  systems  wind  systems,  energy  storage 
systems) 

33  Site  data  and  systems  integration  (energy  resource  data 

energy  consumption  data,  integrating  energy  systems) 

34  ENVIRONMENTAL  PROTECTION 

35  Hazardous  waste  minimization 

36  Restoration  of  installations  (hazardous  waste) 

37  Waste  water  management  and  sanitary  engineering 

38  Oil  pollution  removal  and  recovery 

39  Air  pollution 

44  OCEAN  ENGINEERING 

45  Seafloor  soils  and  foundations 

46  Seafloor  construction  systems  and  operations  (including 

diver  and  manipulator  tools) 

47  Undersea  structures  and  materials 

48  Anchors  and  moorings 

49  Undersea  power  systems,  electromechanical  cables 

and  connectors 

50  Pressure  vessel  facilities 

51  Physical  environment  (including  site  surveying) 

52  Ocean-based  concrete  structures 
54  Undersea  cable  dynamics 

82  NCEL  Guides  &  Abstracts  j~|  None- 
91  Physical  Security  remove  my  name 


